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ABSTRACT 

A high-precision two-dimensional stellar evolution code has been developed 
for studying solar variability due to structural changes produced by varying in- 
ternal magnetic fields of arbitrary configurations. Specifically, we are interested 
in modeling the effects of a dynamo-type field on the detailed internal structure 
and on the global parameters of the Sun. The high precision is required both 
to model very small solar changes (of order of 10~^) and short time scales (or 
order of one year). It is accomplished by using the mass coordinate to replace 
the radial coordinate, by using fixed and adjustable time steps, a realistic stel- 
lar atmosphere, elements diffusion, and by adjusting the grid points. We have 
also built into the code the potential to subsequently include rotation and tur- 
bulence. The current code has been tested for several cases, including its ability 
to reproduce the 1-D results. 

Subject headings: Sun: evolution — Sun: interior — Sun: variability 

1. Introduction 

Modern standard solar models are known to yield the solar structure to an amazing 
degree of precision (see e.g., Guenther & Demarque 1997; Basu, Pinsonneault & Bahcall 
2000; Winnick et al. 2002). These models, however, cannot explain the solar cycle, and 
other solar-cycle related variability. The reason for this shortcoming is that these models do 
not include the dynamo magnetic fields and relevant temporal variability. 

Following the suggestion by Sofia et al. (1979) that any change in the solar lumi- 
nosity, L, must be accompanied by a change in the radius, R, a number of theoretical 
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investigations have attempted to establish the relationship between these changes (denoted 
as W = A In i?/ A In L), by including internal processes designed to mimic the effects of 
dynamo fields. We classify them into three broad categories: 

(i) Perturbation calculation. See Endal et al (1985) for a review of the early work, and 

Balmforth et al. (1996) for subsequent work. 

(ii) Approximation analysis. See Spruit (1991, 2000) for reference. 

(iii) Stellar evolution with magnetic fields. This method was initiated by Lydon & Sofia 
(1995), updated by Li & Sofia (2001), generalized to include turbulence by Li et al. 
(2002), and further generalized to include the interaction between turbulence and mag- 
netic fields by Li et al. (2003). 

The first two are illustrative, but not conclusive. The third can model the effects of 
arbitrary magnetic field configurations. Li et al. (2003) attempted to produce the observed 
cycle variations of all global solar parameters and the p-mode oscillation frequencies. The re- 
sult is promising (e.g., Sofia et al. 2005), but it is not final both because the one-dimensional 
approximation is utilized, and because not all global parameter data exist for the same time 
span. The 1-D approximation only allows us to use a shell-hke magnetic field configura- 
tion. This approximation is relatively limiting. For example, in one-dimensional codes the 
energy flux can only advance to the surface by penetrating the magnetic fleld shell. If the 
magnetic fleld were toroidal, like most dynamo models require, energy flow could circumvent 
the fleld. The aim of this paper is to describe a mathematical technique that can model 
arbitrary magnetic field configurations by generalizing our one-dimensional technique into 
the two-dimensional case. 

In order to match the observed variations of solar global parameters and helioseismic 
frequencies, two-dimensional solar models should fulfill at least the following precision re- 
quirements: 

(a) A luminosity resolution equal or better than 10^^% per year because the observed cyclic 

variation of total solar irradiance is about 0.1% per cycle, 

(b) A radius resolution equal or better than 10~^% per year because the observed cyclic 
variation of solar radius may be as small as 10~^% per cycle, 

(c) A realistic atmosphere model because the helioseismic frequencies are sensitive to it, 

(d) Suitable boundary conditions because the model is sensitive to them. 
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(e) Elements diffusion because the helioseismic frequencies are sensitive to composition, 

(f ) A magnetic field because there is no cyclic variation without magnetic field, 

(g) Turbulence because helioseismic observations require it, 

(h) The interaction between turbulence and magnetic fields because hehoseismic observa- 
tions require it. 

Our one-dimensional code, which is based on the Yale Stellar Evolution Code YREC (Guen- 
ther et al. 1992), meets all these requirements, a nontrivial accomplishment. It is difficult to 
modify the other existing two- or three-dimensional codes (e.g., Dcupree 1990; Turcotte et 
al. 2001) since each of them was developed with specific objectives not requiring this degree 
of precision. 

We attempted to include magnetic fields in Deupree's two-dimensional stellar evolution 
code (1990), but we were unable to compare the model results with solar observations and 
our one-dimensional results, probably because 

(1) the two-dimensional model has different center and surface boundary conditions than 

the one-dimensional model, 

(2) the two-dimensional model does not include an atmosphere model, 

(3) the numerical accuracy is not high enough to match the solar observations. 

This experience convinced us that it would be easier to develop a high precision two- 
dimensional stellar structure and evolution code by straightforward generalizing our one- 
dimensional code rather than modifying an existing two-dimensional code. Our experience 
shows that this conviction was well founded. 

The highest precision requirement is that the cyclic variation of solar radius should 
be better than 10~^% per year because the observed cychc variation of solar radius may 
be as small as 10~^% per cycle. There are various uncertainties in the input physics (e.g., 
Boothroyd & Sackmann 2003; Sackmann & Boothroyd 2003). Though these uncertainties 
affect the interior structure of the Sun, they have little influence on the cyclic variations of 
solar global parameters such as solar radius, solar luminosity, and solar effective temperature 
because of calibration and subtraction of the same parameter at two different times, which 
remove various possible uncertainties in the cychc variations of global solar parameters. Such 
a high precision for the cyclic variations of global solar parameters is thus achievable. 
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We outline here the basic schematic of the method in order to prevent the readers from 
getting lost in the detailed derivations. 

As it is a common practice, the starting points are the conservation laws of mass, mo- 
mentum, energy, and composition, as well as the Newtonian universal gravitational law. 
Both momentum conservation equations and the Poisson equation are second-order differen- 
tial equations. We use the radiation transport equation to relate the temperature gradient 
to the energy flux in the radiative zone, and use the mixing length theory to calculate the 
temperature gradient in the convective zone. We include magnetic fields in this paper, and 
include in the code the potential to subsequently include turbulence and rotation. 

The main relation is the coordinate transformation from the radial coordinate r to the 
mass coordinate m. Regarding mass, we should specify the spatial range that the mass 
occupies. We use the equipotential surface on which 

$(r,^;i) = $e (1) 

to indicate the spatial range, where we have assumed that the system is azimuthally- 
symmetric or axi-symmetric and that <l>c may vary with time. The time coordinate t is 
taken as a parameter. Solving Eq. (1) for r, we obtain the equipotential surface: 

r^R{^c,0;t). (2) 

This equipotential surface encloses volume V$, which is defined by 

r G [0, 2tt] 

V^ = l ee[0,n] (3) 

[ re [o,i?($e,^;t)]. 

The mass contained in V$ is defined by 

pdV = / d(j) d6sm6 / dsp{s,9;t)s^ 
Jo Jo Jo 

= 27r / desine / dsp{s,e-t)s^ , (4) 



where p = p{r, 9; t) is the density. Solving Eq. (4) for $c, we obtain 

$e=*cKi)- (5) 

Substituting Eq. (5) into (2), we obtain the coordinate transformation relation from (r, 9; t) 
to (m, 9; t): 

r^R{^c{m;t),9;t)^r{m,9;t), 9^9, t^t. (6) 
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For any dependent variable X, for example, P, T, Fr, or p, we have 



/dX\ ^ fdX\ fdX\ fdr_^ 



In order to achieve a high precision that is comparable to the one-dimensional solar 
model in the two-dimensional case, using limited computational resources, we cannot directly 
numerically solve those conservation equations and the Poisson equation. For example, even 
in the hydrostatic case, we have five dependent variables such as pressure (P), temperature 
(T), radius (r), gravitational potential ($), and flux {F^ or L = Anr'^Fr). The coeflicient 
matrix of the linearized difference equations with grids M x N has A/" = 5MN x 5MN 
elements, where M (N) is the number of grid points for the mass (co-latitude) coordinate. 
The one-dimensional solar model has M > 2000. If we take N = 20, we obtain J\f > Ax 10^°. 
Since 2^^ ^ 4 x 1024^, a 32-bit computer can handle only 2 x 1024^ ~ 2 x 10^ elements, 
noting that 1 bit is used to represent the sign of a number. Of course, a 64-bit computer 
does not impose such constraint, but the computation speed will become an obstacle. 

Analytical solutions are accurate, but such solutions are hard to obtain in the general 
case. The 1-D case is accurate because we do not need to numerically solve the second-order 
Poisson equation for the gravitational potential $o- It is well-known that the gravitational 
acceleration in the spherically symmetric case is 

g = d^o/dr = Gm/r^. (8) 

In order to take similar advantage in the 2-D case, we show in the paper that Eq. (8) can be 
generahzed as follows 

Gm ^ ^ , . cot ^9$ 

where 0(2) represents the much smaller correction than the retained, and pm is defined by 

p^{m,e-t) = ^^ I deR\^c,0;t)p{R{^c,0;t),e;t)sine. (10) 
Jo 

Like Eq. (8) in 1-D case, Eq. (9) substantially simplifies the 2-D stellar structure equations. 

In the two-dimensional case, the radial component of the energy fiux vector F, Fr, and 
the ^-dependent luminosity, L = 47rr^Fr(r, 9; t), are equivalent to each other, but the actual 
luminosity L* is different than the ^-dependent luminosity L because 

L* = 27r/ Fr{r, 9; t) sin ed9. (11) 
Jo 
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The basic equations are described in section 2, and then the coordinate transformation 
from the radial coordinate to the mass coordinate is performed in section 3. Various possible 
magnetic field configurations are converted into suitable expressions that appear in the stellar 
structure equations in section 4. Boundary conditions are equally important. So we use a 
whole section (section 5) to elaborate them. The method of solution is detailed in section 6. 
The coefficient matrix and input physics used in section 6 are presented in Appendix A and 
B, respectively. The evolution sequences without any magnetic field and with a shell-like 
magnetic field are presented in section 7 and section 8 to test the method. 



2. Basic equations 

The basic equations consist of the timc-dcpcndcnt conservation laws of mass, momen- 
tum, energy, and composition, the Poisson equation (Deupree 1990), as well as the radiative 
transfer equation (Unno & Spiegel 1966): 

§^ + V • (pv) = 0, (12a) 

J 1 

p— = -VP-pV$ + — (V X B) X B, (12b) 
at Att 

pT^ = pe-V-F,,d, (12c) 

^ = Qi^ (12d) 
= AnGp, (12e) 
V-Frad = -4«:p(>^-^), (12f) 

where v is the velocity of a fluid element, B the magnetic field, e the nuclear energy gen- 
eration rate per unit mass, Fj-ad the radiative energy flux, pi the density of species i, Qi the 
creation rate of species i, G the universal gravitational constant, J the mean radiative inten- 
sity, K the absorption coefficient, and B the Kirchhoff-Planck function. The total derivative 
is defined by d/dt = d/dt + v ■ V. 

The specific entropy St includes both nonmagnetic and magnetic components, as shown 
in the first law of thermodynamics (Callen 1996; Lydon & Sofia 1995): 

TdSr ^dU + PdV - dx, (13) 

where U is the nonmagnetic specific internal energy, V = 1/p is the specific volume, 
X — |Bp/87rp is the specific magnetic energy, P is the nonmagnetic pressure. Since the 
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magnetic work dx is taken from the nonmagnetic internal energy, the total internal Ut 
energy decreases: 

Ut^U-x- (14) 

The isotropic magnetic pressure component Pm can be expressed by x P' 

Pm = XP- (15) 

The total isotropic pressure component Pt can thus be defined as follows: 

Pr^P + Pm. (16) 

Using Pt, T and x ^ independent thermodynamic variables, the equation of state and the 
first law of thermodynamics read (Lydon & Sofia 1995): 

dp/p = adPr/ Pt — 5dT /T ~ udx/Xi (17a) 

TdSr = CpdT-{6/p)dPT + {PTSiy/Pma)dx, (17b) 

where 

a = id\np/d\nPT)T,^,t, 5 = -{din p/d\nT)p^,^.„ (18a) 

u = {din p/ din x)Pr,T;u Cp = {OUt /dT)p^^^^. (18b) 

From the first law of thermodynamics (Eq. 17b), we can define two adiabatic gradients, one 
fixes the specific magnetic energy: 

„ f dlnT\ PtS 
another does not fix the specific magnetic energy: 

^^"(^^) =Vad(l-^V», (20) 
\ ^ / St 

where the magnetic energy gradient are defined as follows 

„ 9 In Y , , 

In order to close the radiative transfer equation (Eq. 12f), we use the Eddington ap- 
proximation (Unno & Spiegel 1966): 

Frad = -T^VJ. (22) 
3kp 

Unhke Deupree (1990), we do not directly solve these equations. We first perform some 
analytic work to make some approximations in advance. 
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2.1. Mass conservation equation 

Dcuprce (1990) uses the constancy of the total mass during the model evolution to 
determine the radius at the equator. In contrast, we want to determine the equipotential 
surface S^: r = i?($c,^;t) = r{m,6]t), as in the one-dimensional case. 

Mass conservation can be expressed by either Eq. (4) or its differential form: 

^ = ^ = 47rr2(m, 9; t)p^{m, 9; t), (23) 

where 

r\m, 9; t)pm{m, 9;t) = ^J^ d9R\^,, 9; t)p{R{^„ 9; t), 9; t) sin 9 = f{m; t). (24) 

It should be pointed out that in general, 

Pm{m,9;t)j^p{R{^c:9;t),9;t). (25) 

Nevertheless, in the spherically-symmetric case, Pm{fn; t) is indeed equal to p(i?($c); 0- Since 
f{m;t) is an integral, the two-dimensional case is much more complicated (i.e., nonlocal) 
than its one-dimensional counterpart (local). This complexity may be the price we have to 
pay from one dimension to two dimensions. 



2.2. Gravitational acceleration 

We want to show here that the last two terms [excluding 0(2)] on the right hand side 
of Eq. (9) are due to the 2-D corrections to the gravitational acceleration. To this end we 
should start from the Poisson equation, Eq. (12e), which can be expanded as follows in the 
spherical polar coordinate system: 

where we have assumed that $ = $(r, 9] t) does not vary with the coordinate. We expand 
$ around its spherically-symmetric state: 

$(r,^;i) = $o(r;i) + 5$(r,^;i), (27) 

where 5$ is a small correction, and 
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Substituting Eqs. (27-28) into Eq. (26), we obtain 



Gm ^ ^ . . cot ^9$ 

+ 27rGr{p - - —-— + 0(2), (29) 



dr r2 ^"^^ 2r 09 

where 



2.3. Momentum conservation equation 

Generally, we can decompose the total velocity v in the basic equations into three 
components: 

V = Vo + Vrot + v', (31) 

where Vq is a secular evolution velocity, Vrot the rotation velocity, and v' the turbulent 
convection velocity. We neglect the secular expansion and rotation velocity components in 
the momentum conservation, i.e., we assume 

V = v' (32) 

in Eq. (12b). We checked in the one-dimensional case that the term dVo/dt in the momentum 
equation is negligible. Substituting Eq. (32) into (12b) and averaging the resultant equation 
over the time t and azimuthal angle 0, we obtain 

pVt7^=-VP-pV$ + -^(VxB) xB, (33) 

47r 

where v'"^ = v'J^+v'y'^+v'J^ is computed by solving the basic equations in the three-dimensional 
convection simulations of the outer layers of the Sun (Robinson et al. 2003), in which the 
average is taken over the time t and the horizontal coordinates x and y in a sample box. 
We have shown how to include turbulence in the one-dimensional case (Li et al. 2002). 
We neglect the turbulent contribution to the momentum equation here so as to stress the 
two-dimensional effects due to magnetic fields, i.e., we simply set 

v^^O (34) 

in this paper. 

We assume that the system is azimuthally symmetric. Under this assumption, the vector 
equation (33) is equivalent to the following two scalar equations: 

OPt (9$ , , 

^ = + ^^^^^ 

= — ^ + ^«' 35b) 
r o9 r o9 
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where Pt — P + Pm^s the total pressure, including the magnetic pressure Pm — B^/Sn, and 

7^=i_(B.V)B. (36) 

Noticing that 

1 /B^\ 1 

— (VxB)xB = -V — +— (B-V)B. (37) 



47r VS"^/ 47r 

In the one-dimensional case, we have only a single scalar equation to describe the mo- 
mentum conservation, i.e., Eq. (35a). In contrast, we need three scalar equations for the 
momentum conservation in the two-dimensional case, i.e., Eqs. (29), (35a) and (35b). It 
would be much better if we can combine these three equations into a single scalar equation. 
Fortunately, we can. For this end, solving Eq. (35b) for d^/d9, we obtain 

ad p o9 p 

Then substituting this into Eq. (29), we obtain 

(9$ Gm ^ , , cote OPt cote^, ^, ^ 

_ = _ + 2.Gr(p - p„) + __-_„. + 0(2). (39) 

Substituting Eq. (39) into Eq. (35a), we finally obtain 

^ = + «^ - - p„) - + \h, cot « + 0(2). (40) 

This is our momentum conservation equation. The last three r.h.s. terms represent the 
two-dimensional effects. 



2.4. Energy conservation equation 

The energy conservation equation (Eq. 12c) depends on the velocity in the total deriva- 
tive: 

+ (41) 
The secular expansion velocity Vq cannot be neglected, and from now on we define 

f.^f + V..V... (42) 

The statistical average of pTV ■ VS't, namely < pTv' • V St >, will determine the divergence 
of the convective flux Fconv^ 

V • Feonv =< pTV ■VSt>. (43) 
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By defining the total energy flux to be the sum of both the convective and radiative flux, 
F = Frad + Fconv, Eq. (12c) becomes 



V-F = p e-T 



'dt 



where 



dt 



dlnT 
dt 



ad 



a 



dt 



(44) 



(45) 



In the azimuthal case, Eq. (44) is equivalent to the following equation 



1 djr^Fr 



I r^dS , 



1 djsmeFg) 
r sin 9 dO 



(46) 



We will work out both the radial flux component F^. and polar flux component Fq in the 
next subsections. 



2.5. Energy transport by radiation 

The radiative flux is given by Eq. (22), in which the mean radiative intensity J is 
governed by the radiative transfer equation (Eq. 12f). The Plank function B is known: 

B = ^^T\ (47) 

where a is the radiative constant, and c the speed of light in vacuum. In stellar interior, 
local thermodynamic equilibrium is a good approximation, which leads to 

j^B^ —T\ (48) 
47r 

The more accurate solution of Eqs. (12f) and (22) is (see Unno & Spiegel 1966): 

J^B + i^^B + ^V^B + ■■■, (49) 
3 5 

where Ip = 1/ Kp is the mean free path of photons. Since 



^ r"^ dr \ dr ) ^ ^ ^ 

using Eq. (47) in Eq. (50), we obtain 

V^B^4V.(4V,-l+«-W. + §;^-i^:l^)^, (51) 
\ c^IuPt m J Hp 
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where Hp = —dr/dlnPT — P/ pg is the pressure scale height and V5 is the actual temper- 
ature gradient. Substituting Eq. (51) into Eq. (49), we obtain the mean radiative intensity 
that goes beyond the local thermodynamic equilibrium approximation with one more term 
correction: 



J = 1 + Ao 



/2 
rip 



B, 



where 



An 



4V^ -l + a-5Vs 



d In Pt 



(52) 



(53) 



We want to note that the term {AT^r'^Hppm/m){lp/ H%) ^ 1 is neghgible in the whole star. 
Using this solution in Eq. (22), we obtain 



where 



A= An 



^ 2 



(91n/t 
dhxT 



Pt 



rad 



1 1 

2v;; 



AacT^ 

2)Kp 



[I + A)vr, 



(54) 



1 + 



(91n/t 
dlnPr 



2H 



p 



ld\n\o\ I 



4 dhiT 



Since Ip is much smaller than Hp in the optically-thick region, we know A ~ so that 
Eq. (54) reduces to the widely-used approximation expression without A. However, Ip can 
be comparable to or larger than Hp near the surface, the correction factor A cannot be 
neglected. 



2.6. Energy transport by convection 

Without solving the turbulent convection problem, Eq. (43) only tells us that the con- 
vective flux may depend on the convective velocity Vconv and the entropy St, where the 
convective velocity Vconv has only the statistical meaning. We use the mixing-length theory 
to obtain an analytic expression for Fconv in terms of Vconv and St (e.g., Stix 1989; Lydon & 
Sofia 1995). Since the convective velocity has only the statistical meaning we assume that 
the turbulent convection is isotropic so that Fconv depends on the amplitude of the convective 
velocity, Vconv- 

Fconv = -:^pTlmf{Vconv)^ST, (56) 

where f{v) will be determined by the mixing-length theory, 1^ is the mixing length. It is 
well-known that f{v) = v when the radiative loss of the convective element and the magnetic 
fields are neglected (e.g., Stix 1989). 
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The starting point of the mixing-length theory (MLT) is to calculate the excess heat 
flux in the radial direction: 

-^conv = PVcorwDQ = pVconviQe- Qs) 

= pVconACp{Te - T,) - {6/p){PTe " Pts) + {PtSu / Pr^a){Xe " Xs)], (57) 

where we have used the first law of thermodynamics DQ = TDSt- The subscriptions "c" and 
"s" stand for a convective eddy and its surroundings. If the eddy is always assumed to be in 
pressure equilibrium {DPj- = Pxe = Pts = 0) and magnetic equilibrium {Dx = Xe ~ Xs = 0) 
with its surroundings, we have 

^c^onv = P^convCp(re - T.) = pV, on. CpT{V , - Ve), (58) 

Zll p 

Where the mixing-length approximation in MLT is used to calculate the temperature (or 
density) difference: 

We have also defined the eddy and surrounding temperature gradients and the pressure scale 
height: 

The convective velocity Vcom is generated by the radial buoyancy. The radial buoyancy 
acceleration is related to the density difference: 

^ = -9{Dp/p), (61) 

where g is the gravitational acceleration. For standard MLT, the density difference is related 
to the temperature difference via the equation of state with DPt — and Dx — (see 
Eq. 17a): 

Dp/p = -{DT/T)5 = M(v, - v.). (62) 
We also use the mixing-length approximation to calculate buoyancy acceleration 

where we have assumed that the convective velocity Wconv equals half of the maximum velocity 
((ir/(it)inax- Substituting Eqs. (62-63) into (61), we obtain 

^conv = ^<^(V.-Ve)(08//p). (64) 
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This gives 

^s-^e^ ^<n. (65) 



Substituting this into Eq. (58), we obtain 

F:^^^ = (ApCpT/gU)vl,^. (66) 



Eq. (54) yields 

AacT^ 



Defining a "radiative" gradient 

3KpHpFr 

Vrad - 



we obtain 

AacT^ 



(68) 



- , „ ^ a + A) Vrad. (69) 
We use the energy flux conservation law -^conv+-^rad ~ to constrain the convective velocity: 



ApCpT SnpHp 3 



1 + A glmS AacT^ 



^conv + V. = Vrad. (70) 



2.6.1. Nonmagnetic adiahatic approximation 

When the convective eddy is adiabatic, its temperature gradient equals the adiabatic 
gradient. The nonmagnetic approximation imphes x = 0. Therefore, the temperature 
gradient in a nonmagnetic adiabatic eddy is determined by 

Ve = VU = Vad. (71) 

Eq. (58) thus becomes 

^conv = -^pT^^mW , (72) 

where we have used the following equality: 




Comparing Eq. (72) with the radial component of Eq. (56), we find 

f{v) = V, (74) 
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as stated above. 

Using Eqs. (65) and (71) in Eq. (70), we obtain the cubic equation of the convective 
velocity: 

Vconv + ^TT^conv = Vrad " Vad- (75) 



1 + A glmS 4acT* gl^S 
The convective instabihty condition in the adiabatic approximation is 

Vrad > V, > Ve = Vad, (76) 

according to Eq. (64). 



2.6.2. Nonmagnetic nonadiabatic approximation 

During its rise the eddy radiates energy into its environment. For this reason the eddy 
gradient Ve differs from the adiabatic gradient, Vad- We decompose the convective flux 
(Eq. 58) into the adiabatic (the first r.h.s. term) and nonadiabatic (the second r.h.s term) 
fluxes: 

1 1 V 

i^^onv = ^^pTCp(V.-Ve) 

= ^^^pTCp(V.-Vad) + ^^^pTCp(Vad-Ve). (77) 

If the effective cross section of the convective eddy is q, the heat energy loss rate of the eddy 
due to radiation can be expressed by 

"df ^ 2 Hp P^^P^^^ ~ ^^>^- 
We can also use Eq. (54) to calculate the radiative loss: 

dQr ^ AacT^T,-T ^ 1 Z^S 4acT^ _ 

dt Znp d 2Hpd 3kp ^ * ^ ^ 

where d is the effective radius of the eddy, E is the eddy surface. Comparing Eq. (79) with 
(78), we obtain 

Ve - Vad = (^^o/ W) (V, - Ve) , (80) 

where 

/„S4acT3 1 

^0 = , o -, • (81) 

qd SpCp IjnK'P 
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Substituting Eq. (80) into (77), we can express Vs — Ve by Vs — Vad^ 

V. - Ve = \ iy. - V.d) = ^ (^) , (82) 

1 + t;o/^'conv 1 + Vo/Vconv C p \ dv J ^ 

where we have used Eq. (73). Finally, using Eq. (82) in (77) we obtain 



This shows 



f(v) = (84) 



Using Eq. (65) in (82), we obtain 



8Hp 2 



V, - Vad = ^^^^conv(l + '^o/Vconv)- (85) 

Substituting into Eq. (70), we obtain the cubic equation of the convective velocity: 
ApCpTSKpHp o 8Hp 2 , , , ^ ^ 

-^1^'-- + ^s'--^^ + = - ^^^^ 

The convective instability condition in the nonmagnetic nonadiabatic approximation is 

Vrad > V, > Ve > Vad, (87) 

according to Eq. (64). 



2.6.3. General case 
When magnetic fields are present, we have 




We decompose the convective flux (Eq. 58) into the adiabatic (the first r.h.s. term) and 
nonadiabatic (the second r.h.s term) fluxes: 

i^^onv = ^^PTCP(V.-Ve) 

= ^^prcp(v.-v;,) + ^^prcp(v;,-ve). (89) 
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The heat energy loss rate of the eddy due to radiation now can be expressed by 



dt 2 Hp 



pTCp{V^ - Ve)g. (90) 



The radiation loss rate calculated by Eq. (54) is the same as given in Eq. (79). Comparing 
Eq. (79) with (90), we obtain 

Ve - V',d = ivo/v,ony) (V, - Ve) ■ (91) 

Substituting Eq. (91) into (89), we can express — Ve by Vs — V^^: 

V. - V. = ,^ ' (V. - VL,) = ' ^ (^) , (92) 

1 + Vo/Vconv 1 + Vo/Vconv Cp \ Or J ^ 

where we have used Eq. (88). Finally, substituting Eq. (92) into (89) we obtain 

which leads up to Eq. (84). 

Using Eq. (65) in (92), we obtain 



V. - = ^.Lv 1 + • (94) 



^Hp 2 A , ^0 

Substituting into Eq. (70), we obtain the cubic equation of the convective velocity in a 
magnetic system: 

ApCpT?,KpHp 3 8i7p 2 Vq\ , 



The convective instability condition in the magnetic nonadiabatic case is 

Vrad > Vs > Ve > Vad, (96) 

according to Eq. (64). 

Eq. (95) can be rewritten as follows 

2Aoy^ + Vy^ + V^y-V^0, (97) 

where we have defined the dimensionless variable 

y = VVconv/vo (98) 
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and the following dimensionless parameters 

t'o = 6acT^/pCplmKp, 



c 



8Hp ' 



9 1 



81 + A 

We choose Im^/qd — 9/2 for spherical eddies and d/lm — 8/9. The convective gradient can 
be expressed by y: 

Vconv = V, = + (V,ad - VJy{y + V), (99) 

according to Eq. (94). 

When magnetic fields are neglected, — Vad, all formulas automatically reduce to 
their counterparts in section 2.6.2. 



2.7. Energy flux vector 

In the convective zone, the total energy flux vector equals the radiative flux (Eq. 54): 

F ^ + X)VT. (100) 

Sup 

In the convective zone, the total energy flux vector equals the sum of the radiative (Eq. 54) 
and convective (Eq. 56) fluxes: 

F = -1^(1 + A) Vr-- ^^^"-^^""^ VSt 
Znp 2 1 + t;conv/t'o 



AacT^ (1 _^ A) + - P^P^""^" 



I pCpTVjrnVcon. ^ ^^ ^^ 

2 1 + w/^^o Pt ^ ^ 



Sup 2 1 + t;conv/t'0. 

where we have used the following formula 

VSt = {Cp/T)VT - {CpVJPt)VPt. (102) 

2.8. Composition conservation 

Eq. (12d) describes the composition conservation law, which can be rewritten as follows: 

P^ + X,^ + V- {pXy) = Q„ (103) 



-19- 



where we have used Eq. (32) and have used the mass fraction Xi = pi/p to replace density pi. 
We have also assumed V- v' = 0. Eq. (103) involves two timescales: one is the thermonuclear 
reaction timescale Tnuch which determines Qi and is quite long, another is the convection 
timescale Tconv) which determines the convection mixing and is much shorter than the former. 

As before, taking the statistically average over Eq. (103), we obtain 

+ ^V- < pXy g,, (104) 

where we have used the assumption < dp/dt >= and defined = Qi/ p. Using the mixing 
length theory, we can express the mass flux Fj =< pXjv' > as follows: 

Fj = —pV^nylmSXi. (105) 

Substituting Eq. (105) into Eq. (104), we obtain 

dX 1 

= + ■ (pt^conv/mVXi). (106) 

In the radiative zone, the clement diffusion velocity Wi (e.g., Thoul et al. 1994) changes 
the local composition in addition to the thermonuclear reactions. Element diffusion in stars 
is driven by pressure gradients (or gravity), temperature gradients, composition gradients, 
and radiation pressure. Gravity tends to concentrate the heavier elements toward the center 
of the star. Temperature gradients lead to thermal diffusion, which tends to concentrate 
more highly charged and more massive species toward the hottest region of the star, its 
center. Concentration gradients oppose the above two processes. Radiation pressure causes 
negligible diffusion in the solar core. Element diffusion affects the element abundances, the 
mean molecular weight, and the radiative opacity in the radiative zone, and therefore affects 
the calculated neutrino fluxes and oscillation frequencies, on which observations impose strict 
constraints on the solar model. 

The characteristic time for elements to diffuse a solar radius under solar conditions 
is of the order of 6 x 10^^ yr, much larger than the age of the Sun. Element diffusion 
therefore introduces only a small correction. Many authors have studied this topic carefully 
(see Thoul et al. 1994 and references therein), and both portable subroutine and analytic 
formulae for element diffusion calculations are available. In particular, the formulae for the 
element diffusion velocity flts our theoretical framework developed in this paper. We use the 
formula given by Thoul et al. (1994) with included: 
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where 

See Thoul et al. (1994) for the expansion coefficients, which are actually computed by 
numerically solving the multi-fluid equations for all species. These formulae just give readers 
the main idea. We use the portable subroutine provided by the authors to compute the 
element diffusion correction. Diffusion in the polar direction is negligible. 



3. Coordinate transformation from r to m 

So far, all derivatives with respect to 9 assume r to be constant. What we need is 
to obtain the corresponding derivatives at the constant m. This can be done by using the 
so-call implicit-function rule, that is 

d_\ _ rdr\ d_ _ / d\nr \ d 

d0j^~\d9j^^[d9)^dr-[d9)^^[ 09 J^dlnr ^ ^ 

Prom now on, we use the following shortcuts to save writings: 

r' = In r, p' = In p, P' = In Pt, T' = In T, and s = In m. (110) 

We note that In is the natural logarithm. 

Another formulas we need for this purpose is the mass conservation equation, Eq. (23), 
which can be rewritten as follows: 

dr' m ^^^^^ 



ds Anr^pr, 



3.1. Momentum conservation equation 



We perform the necessary coordinate transformation from r to m in Eq. (40). The only 
term that needs to be transformed is the term that contains dPr/dO, which is equivalent to 
{dPT/d9)r. Using Eq. (109), we obtain 



dP, 



T 



d9 



dPj 
"89 



dp 



dr' \d9 J, 



(112) 



Consequently, Eq. (40) becomes 

dPr 
dr 



+Hr + -He cot 9 



cot 9 


(dr'\ ' 


-1 - 


2 







Gmp 



2TTGrp{p - pm) 



cot 9 fdP, 



2r 



89 



+ 0(2). 



(113) 
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The first factor on r.h.s. is caused by the coordinate transformation from r to m. 



Since dPx/dr = [Pt / r){ds / dr'){dP' / ds) , using Eq. (Ill), we can rewritten Eq. (113) 
as follows: 



where 



e 



M 



dP' 

ds 



Gm? p 



47rr4Pr P, 
Gm{p-pm) p 



+ e + M + 0{2), 



'^rPr Pm 
m cote fdP' 



1 - 



cot^ / dr' 



dcote 



09 



^nr^Pm 2 
Gm'^ p cot 6' / dr' 

At^t^Pt'^'^ \de 



cot 9 f dr'\ 
99), 
cot 9 f dr'\ 



m 



ATcr^PmPr 



1 - 



cot 9 ( dr' 



89 



-1 



{Hr + -He cot ( 



(114) 



3.2. Energy conservation equation 



The starting equation is Eq. (46). The only term that needs to be transformed is the 
term that contains the derivative of (sin^F^t) with respect to 9. This term is a small 2- 
D correction to the Energy conservation equation since Fg, which is given in Eq. (101), is 
already a combination of the first-order derivatives of T and Py, 



AacT^ ,^ IpCpTlmVc 
— 1 + A + ^ ^ 

6Kp 



2 1 + ^^conv/'l^O 



r 89 2 1 + Vconv/vo r 89 ' 



(116) 



Therefore, after neglecting the higher-order corrections as we did above, the energy conser- 
vation equation becomes 



1 8{r^Fr 

J.2 Qj. 



pie-T 



dSj 
~dt 



Fe cot 9 



+ 0(2). 



(117) 



This shows that we only need to transform Fq from r to m. Applying Eq. (109) to (8T' /89)j 
and {8P'/89)r in Eq. (116), we obtain 



Fe 



AacT\ IpCpTlmVc 



2 1 + t'conv/t'O 



dr 

'89 



dr 

8r' 



dr' 
'89 
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+ 



1 pCpTlmVconv ^'ad 



2 1 -|- Vconv 

AacT^ 



dP' 



Sup 



/vq r 
1 + A) + 1'^^^'^^™^= 



Qpf (Q^J 



dr' 
1 



2 1 



+ 



1 pCpTlmVconv ^'a.d 



2 1+ t^conv/^'o r 
where V is the temperature gradient. 



dP' 



/vq \ r l\ do 
Gmp f dr' 



+ 



Gmp f dr' 
rPr \~de 



+ 



rPj 



09 



(118) 



The second step is to use Eq. (Ill) to replace dr by ds in Eq. (46). Unlike r', P', and 
T', which are the natural logarithms, we define 



L' = Anr^Fr/LQ, 
which is not a logarithm at all. The resultant equation is 



dL' 
ds 



^—m e — T^— — 



^0 



dt J p. 



P 



1 mFa cot 6 



^0 



rpv 



+ 0{2). 



(119) 



(120) 



3.3. Composition conservation 

Eq. (106) involves the derivatives with respect to 6 at constant r. Since what we need are 
the corresponding derivatives at constant m, this equation needs a coordinate transformation 
from r to m. To this end, we first expand it as follows: 

1 d f o , dXA 1 d f . ^ , dXA 

+ r P^^conv^m^ + TTT-— 7^ 1^ OpV^^n^lm^ ■ (121) 



ot 2r^p or \ or J Ir'^p sm U oO \ oU 

The last rhs term needs the transformation. We retain its most important part. The resultant 
equation is 

( ldp'\ dXi cot ev^onJm f dXA .^i^. i^^^. 

where 

^/gx_l ^^^» Vconylm dp' dXj VcorJmd'^Xi COt 9Vconvlm dXj Z' dr' \ 

^ >~ 2 dr^ 2r2 5^ dO 2r^ d9^ 2r'^ dr' \d9 J ^' ^ ' 

We have takan advantage of the fact that I'conv^m ~ constant in stars. 

Since the convection timescale is much shorter than the evolution timescale, the con- 
vection zone is well-mixed on the evolution timescale. As a result, the detailed expression 
for the composition conservation equation in the convection zone does not matter much. We 
do it here just for the sake of completeness. 
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3.4. Two-dimensional stellar structure equations 

In summary, we obtain the two-dimensional stellar structure equations with (m, 9) as 
independent variables: 



dr' 
ds 
dP' 

ds 

dr 

ds 

dv 

ds 

dXj 
dt 



m 



P 



Ai^r^p \ p. 



P 



dP' j Vrad radiative 
ds 1 Vr convcctive 



+ {e} + {M} + 0{2), 



1 f rj^dSr 

— m \ e — 1 

Lq V dt 



_P_ 

Prn 



1 mFa cot 9 



Lq rp„ 



0{2), 



(124a) 
(124b) 
(124c) 
(124d) 



1 + Idp^] V Z ^ + i—V 
'2 dr' I conv''m Qr ~ \ 2r'^ "^conv 



,.2p Qy. (^ P-^i^i) 



cot 6» „ 



(radiative) 
{W)J + 0(2)- (convectiieY^''^ 



Those terms or factors associated with 2-D effects are indicated by {} in the equations. The 
symbols used above are defined as follows: 



e 



M 

Fe 



m cot 9 fdP' 



Anr^Pm 2 \d9 
Gm{p-pra) p 



1 - 



cot 9 ( dr' 



2rP, 



T 



Pr, 



cot 9 f dr' 



d9 
-1 



Gm^ p cot 9 (dr'\ 



cot 9 ( dr' 



_ cot9 (dr'\ 



m 



Anr^PmPr 



d9 



{Hr^ -He cot 9), 



(125a) 
(125b) 



?>Kp 



1 + A) + 



1 pCpTl^v, 



^ 1 pCpTlrnVconv'^'ad 

2 l + T;convM r 



2 1 +^;conv/^^0 

dp' 



dr 



d9 



d9 

Gmp / dr' 
^ rPr \d9 



Gmp / dr'\ 



(125c) 



Prom now on, we'll drop the subscript m in the derivatives such as {dr'/d9)m, 



dr' 
'd9 



dr' 
'd9' 



(126) 



Wherever needed, we'll specify the subscript m or r to avoid confusion. 
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4. Magnetic fields 

Our strategy is to take advantage of analytical results as much as possible. For this 
purpose, in this section we work out the explicit expressions for the terms associated with 
magnetic fields. 

Generally, a magnetic field has three components. Using the spherical coordinate sys- 
tem, it can be expressed by 

B= {Br,B0,B^). (127) 

All three components arc functions of m and 9 in the azimuthally-symmetric case treated in 
this paper. The B-related terms are expressed by 7i, see Eq. (36), which can be expanded 
as follows 

inrn = er{rB-VBr- BeBe- B^B^) 

+ee (rB • VBe - B^B^ cot 9 + BeBr) 

+e^ (rB • VB^ + B^B^ - B^Be cot 9) . (128) 



Consequently, we see 



Ai^rHr = rB-WBr-BeBo-B^B^, (129a) 
AnrHe = rB ■ V Be - B^B^ cot 9 + BeBr, (129b) 
AnrH^ = rB ■ WB^ + B^B^ - B^Be cot 9. (129c) 



We use B to define three stellar magnetic parameters, in addition to the conventional 
stellar parameters such as pressure, temperature, radius and luminosity. The first magnetic 
parameter is the magnetic kinetic energy per unit mass, %, 

X = SV(87rp). (130) 

The second is the heat index due to the magnetic field, or the ratio of the magnetic pressure 
in the radial direction to the magnetic energy density, 7 — 1, 

^^1 + {B',+Bl)/B\ (131) 

Lydon and Sofia (1995) introduced the first two in the one-dimensional case. Here we 
introduce the third one, the ratio of the magnetic pressure in the co-latitude direction to the 
magnetic energy density, -d — 1, 



1? = 1 + (S? + B^)/B\ 



(132) 
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We can use these three magnetic parameters to express three components of a magnetic field 
as follows: 

= [87^(2-7)xp]^/^ (133a) 
Be = [87r(2-^)xp]'/', (133b) 
B^ = [87r(7 + i?-3)xp]'/'. (133c) 



We discuss below various possible cases. Note that any case should satisfy the restriction: 

V • B = 0. (134) 



4.1. B = (0,0,0) 

In this case, 

X = 0, ^? = 1, 7 = 1, 7i = 0. (135) 
Consequently, the term associated with magnetic fields vanishes, namely, 

M = 0. (136) 



Defining 



= (137a) 



we can rewrite © as follows 



cot 9 dr' \ ^ 



1 - 



2 de 

cot 9 dr'\ ~^ 



2 d9 



„,cot^9P'/ cot^9r' , , 



Solving this case will provide us with a standard 2-D stellar model. 
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4.2. B = (0,0,5^) 

Since is assumed to depend upon only r and 9, Eq. (134) is satisfied for any arbitrary 
function = B^{r, 0). In this case, since 

X = i^|/(87rp), 

^ = 2, 

7 = 2, 

we have 

Br = 0, (139a) 

Be = 0, (139b) 

= (Sttxp)'/'. (139c) 
Substituting them into Eqs. (129a-129b), we obtain 

Hr = -2xp/r, (140a) 
Hg = //^cot^. (140b) 

Substituting them into Eq. (125b), we obtain 

A.^-B^(l4cot^«,(l-^^)-\ (141) 



where 



Eq. (134) requires 



This leads to 



4 _ m XP .^42) 

4.3. B= (0,5^,0) 
a(sin^5e) 



80 



0. (143) 



= 5(r)/sin^, (144) 

where B(r) is an arbitrary function of r. 
In this case, since 

X = S,V(87rp), 
^ = 1, 
7 = 2, 
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we have 



Be 



0, 
0. 



(145a) 
(145b) 
(145c) 



Eq. (144) requires that B — (SttxpY^^ sm9 does not depend upon 9. 



In order to calculate A4, we have to calculate Hr and He first. Substituting Eqs. (145a- 
145c) into Eqs. (129a-129b), we obtain 



1 1 



47rr 2 
XP 



09 



dp' dx' 
89 89 



ATrr^Pm { 8p' 8x'\ { 8r'\ 



m 



+ 



8s 8s J \ 89 J 



- 



- 



Substituting them into Eq. (125b), we obtain A4: 



M 



coi9 f dp' dx' 
~ ~T~ \ d9 ^ '89 
8p' 8x!\ cot 9 8r' 
'ds^~ds) ir~d9 



1 - 



cot 9 dr' 



2 89 



1 - 



cot 9 8r' 



-1 



where 



B' = xp/Pt. 



(146a) 
(146b) 

(146c) 



(147) 
(148) 



In this case, since 



we have 



4.4. B = {Br, 0, 0) 



X = B^JiSnp), 
^ = 2, 
7 = 1, 



Br = (87rxp)'/^ (149a) 
Be = 0, (149b) 
B^ = 0. (149c) 
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Eq. (134) requires 

= 0- (150) 

This leads to 

Br^B{e)/r', (151) 
where B{9) is an arbitrary function of 9. Therefore, we know 

B = {SnxpY^'r' (152) 

varies with only 9. 

Substituting Eqs. (149a-149c) into Eqs. (129a-129b), we obtain 

^ 47rrV^ XP f dp' ^ dx'\ 
^ m r \ds ds J ^ 

He = 0. 

Substituting them into Eq. (125b), we obtain M.: 



In this case, since 



we have 



4.5. B = (0, Be, B^) 



X = {Bi + BD/i^Trp), 
Bl 

1? = 1 + 



7 = 2, 



Br = 0, (154a) 
Be = [Sn{2 - d)xpY'\ (154b) 
B^ = [STi{d-l)xpY'\ (154c) 
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Eq. (134) requires 

^^0. (155) 

This leads to 

= 5(r)/sin^, (156) 

where B{r) is an arbitrary function of r. Therefore, we have the constraint that [87r(2 — 
'&)xpV^'^ sin 6' depend only upon r. 

Substituting Eqs. (154a-154c) into Eqs. (129a-129b), we obtain 

Hr = -2xp/r, (157a) 
He = -^(^-l)cot^ + l|^[(2-^)Xp] 

(i? - 1) cot ^ + i I ^ [(2 - i^)xp\ 



r r \ 89 

1 Airr^p^ f dr'\ d 



r m \d9 I ds 
Substituting these expressions into Eq. (125b), we obtain 

,g f ^ cot 9 dr' 



2-^)xp]. (157b) 



M = -B^ 1- 



2 d9 



„^cot9 fdx' dp' dr\ ( cot9dr' 



where 



89 89 89 J \ 2 89 
.8x' ^ 8p' 8r\ cot98r' ( cot98r'\ 

^^^^s^^)^mV-^m^ ' ^^^^^ 



-B^ = B% + -{^-l)cot'9], 
2 



B^ = B\2-'&), 

-87 - 

^og(2-.). 



89 89 
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In this case, since 



4.6. B^{Br,B0,O) 



^ = 1+ 



7 = 1 + 



Bl 



Bl + B, 



2 ' 



we have 



We have used the fact that 



B, = [87^(2-7)xp]^/^ (159a) 
Be = [87r(2-^)xp]^/^ (159b) 
B^ = 0. (159c) 



7 + 1? = 3. (160) 



A meaningful magnetic field should satisfy Eq. (134). For example, r^B^ does not vary 
with r and smOBe does not vary with 0. 

Substituting Eqs. (159a-159c) into Eqs. (129a-129b), we obtain 

^ _(2-^)^ + (2-7)^^^f^ + ??^+^^" 



m r \ ds ds ds 



1/2 XP [ dx' , dp' ^ di' 



m 86 \ds ds ds 

He = [(2-7)(2-^)]V^^ 

r 

+ [(2 - 7)(2 - + ¥ + 



m T \ds ds ds 

^■.XP f dy' dp' dd" 

+p-tf)^(^ + ^ + ^ 

-(0 _ XP ^T^r^Pm dr' f^.^. 9^" 



r m do \ds ds ds 
Substituting these expressions into Eq. (125b), we obtain 



2 de \ds ds ds \ 2 d9 
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where 



de \ds ds ds J \ 2 de 



de de de J \ 2 de 

+513/'^ + ^ , f. cot edr' 



ds ds ds J \ 2 de 



ds ds ds J 2 de \ 2 de 

igcot^ fdx' , , cote dr'Y^ 



-1 



= Bn(2-^)-^[(2-7)(2-^)]n 

gio ^ ^35^2-7), 

iJii = i35[(2-^9)(2-7)]^/2, 

B'' = ^B^[(2-^)(2-7)]^/^ 

^ l-B^cote[{2-'&){2-j)]^/^, 

B^^ = B\ 

97" a , , 



In this case, since 



4.7. B = (£„0,£^) 



X = (5,^ + i?,^)/(87rp), 

^ = 2, 

- 1 ^ ^1 

^ Bl + BV 



we have 



= [87^(2-7)xp]^/^ 



(163a) 
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Bg = 0, (163b) 
= [87r(7-l)xp]'/'. (163c) 

A meaningful magnetic field should satisfy Eq. (134), which requires that r^Br not to vary 
with r. 

Substituting Eqs. (163a-163c) into Eqs. (129a-129b), wc obtain 

^ _(^_i)^+i![!f^^(2-^)f|^+f:+^^" 



r m r \ds ds ds 

He = -(7-l)^cot^. 

r 

Substituting these expressions into Eq. (125b), we obtain 

^._,.(,.--)-V,.(^,|:,^) (,_--)-, 

where 



4.8. B = [Br-, Bq, B^) 

This is the general case, in which all magnetic field parameters x, d and 7 are variables. 
Therefore, we use the general expressions for i?^, Bq and B^ given at the beginning of this 
section. Substituting Eqs. (133a-133c) into Eqs. (129a-129b) to calculate and Hq, we 
obtain 



m r \ ds ds ds 



+[(2-7)(2-^)]V2^|'^ + ^ + ^'l 



y^iirr^PmXPdr' f dx' , dp' , ^7" 



H) = [(2-7)(2-tf)]''25^-(7 + tf-3)?^cot# 



+[(2-T)(2-*)l'/^i^»^f|^ 

m r \ ds 

n^XP f dp' dx' dd" 

+p-<(^ + ^ + i^ 

-10- q^ 47rr XP dr' f dp' dx' d^" 
™ „ az) a„ + a„ + 



9s 9s 



m r d9 \ ds ds ds 
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Substituting these expressions into Eq. (125b), we obtain 



2 d9 J \ds ds ds J \ 2 d9 



89 \ds ds ds J \ 2 d9 



09 09 09 J \ 2 09 



Os Os Os J \ 2 09 
i'Oy' Op' Oi!}"\ cot 9 Or' ( cot 9 Or' 



9s as 9s 7 2 V 2 9^ 

15 cot ^ fOx' , <9p' , dtr\ cote Or' 



-1 



where 



= B^{7 - 1 - ^ cot ^[(2 - 7)(2 - ^)]^/^ + (7 + 1? - 3)^ cot' (166) 



Reahstic magnetic fields in the stellar interior should satisfy the Maxwell equations. 
One of them is the divergence-free condition specified by Eq. (134). Using the coordinate 
(m, 9), this equation reads 

Anr pm sin 9 0{r'^ Br) ^ O{sm9B0) ^ ^ 
m Os 09 

Assuming — C{m) cos^/r', we obtain by solving this equation for BQ{r, 9), 

D / a\ 4.7rrpmdC{s) . 

Beir, 9) — sm 9. (168) 

m as 

So far, we have finished the coordinate transformation from (r, 6) to (m, 9). This allows 
us to use the analytical formula, for instance, and A4, to describe the two-dimensional 
effects. This effort has at least the following rewards 

a) We can control the approximations by neglecting certain terms. 

b) We can understand if certain factor(s) play an important role by including or excluding 

the corresponding term(s) in the numerical calculations. 

c) We may use the existing technique to numerically solve the two-dimensional stellar 

structure equations. 
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d) We can use the analytical expressions to calculate the matrix element coefficients for the 
linearization correction equations. 

We'll make use of these advantages below. 



5. Boundary conditions 

As usual in mathematical physics, the boundary conditions constitute a serious part 
of the whole problem, and their influence on the solutions is not easy to foresee. In the 
one-dimensional stellar model calculations, the boundary conditions cannot be specified at 
one end of the interval [0, Mtot] only, but rather are split into some that are given at the 
center and some near the surface of the star. The central conditions are simple, whereas the 
surface conditions involve observable quantities. The boundaries in the angular direction are 
located at ^ = and either 6 = 7r/2 or ^ = tt. We'll follow Deupree (1990) in using symmetry 
conditions to determine them. Otherwise, the treatment of the boundary conditions is as 
described in Prather (1976) and as implemented in YREC (Pinsonneault 1988). 



5.1. Central conditions 

Two boundary conditions can be specified for the center, defined by m = 0: 

m = 0: r = 0, L = 0. (169) 

Rewriting Eq. (124a) as follows 

3 

dr^ = dm, (170) 

we can integrate it over a small mass interval [0, m] in which = can be considered to 
be constant. The result 

r-f-^V'm^/^ (171) 

V47rp„c/ 

can be considered to be the first term in a series expansion of r around m = 0. Taking the 
logarithm, we obtain 

r' = i[5-log(4W3)]. (172) 
A corresponding integration of Eq. (124d) yields 

L'^^ie- T^) -P- - ^^cote, (173) 
Lq \ dt J pm LQrpm 
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In both cases we have used the proper boundary conditions (169) by taking the lower hmit 
of integration to be zero. 

Eqs. (172-173) are two central boundary conditions that are equivalent to Eq. (169). 

5.2. Surface boundary conditions 

5.2.1. 1-D surface boundary conditions 

Nothing is a priori known about the central values of pressure Pc and temperature Tc- 
So we need to define the surface and specify the surface values of pressure and temperature. 

In principle, we can use a definition for the surface such as 

m = Mtot. (174) 

However, since near the surface m does not change much, this definition is not accurate 
enough. The theory of stellar atmospheres suggests the use of the photosphere, from where 
the bulk of the radiation is emitted into space: 

T = Teff, (175) 

where T^s is the effective temperature. The optical depth of the overlying layers, 

poo 

T — I Kpdr 
Jr 

is equal to 2/3 for the Eddington approximation, 

T'-\T^^{r+^^. (176) 

where R is the total stellar radius. In contrast, the optical depth = 0.312155 of the 
overlying layers is different than 2/3 if the atmosphere is assumed to obey a scaled solar 
T{t) relation given by Krishna-Swamy (1966) 

T^(r) = ^7;^jf[T + 1.39 - 0.815 exp(-2.54T) - 0.025 exp(-30.0T)]. (177) 

Since Teg is the temperature of that black body that yields the same surface flux of energy 
as the star, then 

m = Mtot : = ATrR^aT^s: (178) 
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where a — ac/A: is the Stefan- Boltzmann constant of radiation, Ls is the total luminosity. 
This is one of two surface boundary conditions. 

The second surface boundary condition is the hydrostatic equilibrium condition: the 
pressure at the surface is given by the weight of the matter above. We can well approximate 
the gravitational acceleration by the constant value qq — GMtot/F^, since the bulk of the 
matter above the surface is anyway very close to the photosphere. We hence have 



is calculated in the following way: the starting values of (Pq) '^o) are chosen by selecting a 
small density po and then computing, 



where Tq = T{t = 0). Then (Pq, Tq) gives pi which gives k,o{pi,Tq) which gives ri = KoPq/q 
or 5t — Ti — Tq. Thus we have Iq — 8t/kq. Then redefining tq = Ti and po — Pi- This 
method could be iterated upon by redefining Tq — T{tq) and so forth: 



Sufficient accuracy was achieved in the atmosphere integration by choosing a small enough 
po (e.g., po = 10"^°) such that r < 10"^. 

From the calculation description, we can see that / = I{Ps,Tcs). However, we do not 
know the explicit expression. Therefore, we cannot directly use Eqs. (178) and (179) as our 
surface boundary conditions. Instead, we solve the following system (Kippenhahn 1967), 




(179) 



where the integration 




Po = (a/3)ro' + Po^ro, 





(180) 



for the a[s which are used for the surface boundary conditions. 



E! = aiP' + aaT' + as, 
InL = a4P' + agT' + ae. 



(181a) 
(181b) 



and for the calculation of the effective temperature. 



= ayP' + agT' + ag. 



(182) 
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Here, the {P', T') refer to the values at the outermost mass point in the model. The last three 
equations can be considered to be the first term in the series expansions of Eqs. (178-179). 

The initial model with an estimated (InL*, T'^q) is triangulated in the (InL, Tgg)-plane 
by constructing three atmospheres of the form 

Al: (lnL*-|Ai,T;|+lAr) 
A2: (lnL*-iAi,r;^-iAT) 
A3: (lnL* + |Ai,r;^). 

If subsequent models or if the model itself during convergence moves significantly out of the 
triangle, the triangle is fiipped until it once again constrains the model. The decision as to 
which point of the triangle should be fiipped (if any) can be made by testing, 

c, = /[(lnL,+i - lnL,+2)(T,'ff - U^^J + (T^'^^^, - T^'^,^^ J (In L - lnL,+i)], 

where / = ±1 is the orientation of the triangle (e.g., in the example given / = +1) and 
+ l,i + 2} is {123}, {231} or {312}. The value of q is tested against eA/^A^ where 
setting e = gives exact triangulation and e > allows the point (In L,Tgg) to be at most e 
of a triangle outside. Begin testing with i = 1 to 3, if c, < — eA^A^^ then fiip point i. 

In Li ^ In Lj+i + In - In Lj, 

rjif J rpl I rp! rjif 

/ ^ -/ 

and repeat the testing again starting with i = 1 until Cj passes for i = 1 to 3. The atmospheres 
that have been fiipped are then recomputed as are all the coefficients Oj. 

This treatment of the surface boundary conditions is the same as that in one-dimensional 
model calculation, except that we move the fitting point to the surface where T — T^s. There- 
fore, we do not need an envelope integration. This has been tested for the one-dimensional 
model calculations and it turns out to be acceptable. This saves much computation time in 
the two-dimensional case. 

Our surface boundary conditions are much more complicated than Deupree's (1990) 
because our applications to the Sun are very sensitive to the surface conditions. 

5.2.2. Deupree's 2-D surface boundary conditions (Deupree 1990) 

In his two-dimensional rotational models, Deupree (1990) uses the following surface 
boundary conditions: 

P — Pref, T — Tref, 
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where pref and Tref are reference density and temperature, respectively. The most difficult 
part of using these surface boundary condition is how to select the reference density and 
temperature at the surface. 

Unlike Deupree, we use Pt and T as independent thermodynamical variables. Since 
p — p{Pt,T), the equivalent surface boundary condition is 



ref ) 



T = T,ef. (183) 



In order to compare this with the standard surface boundary condition given above, we use 
the surface values of Pt and T obtained by using the standard surface boundary condition 
for the current Sun as the reference. Fig. 1 shows the reference values as functions of age 
and their polynomial fits. The fitting formulas are: 



p _ 1 73695.514- 9004.5498-t+ 13898.511-^2 < t < 0.27 , , 

ref - i ^2777.060 - 2211.7088 ■ t - 49.075155 • 0.27 < t < 4.55 ^ 



^ _ . 5647.8836 + 266.07365 ■ t - 539.35360 ■ < t < 0.27 , . 

~ ' 5673.6126 + 28.625469 -t- 1.1516435-^2 0.27 < t < 4.55 ^ ^ 



The age t is in giga year. 



5.3. Polar boundary conditions 

Eqs. (124b-124e) are singular at the poles = and tt) because of ©, A4, Tq, and 
U0. However, if dr' /dO — 0, the singularity due to © will disappear. In order to guarantee 
dr' / do — 0, we also need to zero the other derivatives. Therefore, we require 

dr' _dL' _dP' _dT' _ 

at the poles. In order to remove the singularity due to Al, we have to zero x at the poles, 
namely, 

X = (186) 

at the poles. Eqs. (185-186) arc the polar boundary conditions. Eq. (185) are similar to 
Deupree's (1990) polar boundary conditions, which are the symmetry conditions. 



5.4. Equatorial conditions 



Eqs. (124a-124d) show that the two-dimensional stellar structure equations are not 
singular at the equator. Therefore, there are no special constraints there. If we neglect 



-39- 



0{2) in Eqs. (124b, 124c), the two-dimensional stellar structure equations are a set of first- 
order differential equations. Since we have already specified four boundary conditions at 
the north pole {9 = 0), wc do not need extra boundary conditions at the equator. If we 
want to include those terms that contain those second-order derivatives in 0(2), we have to 
specify four equatorial boundary conditions or five polar boundary conditions at the south 
pole (6* = tt). We do not include those second-order derivatives in 0(2) in this paper for the 
following reasons: 

1) They are much smaller corrections than the retained; 

2) They may cause a much bigger numerical error than the actual corrections; 

3) They require a totally different method of solution (e.g., Deupree 1990). 



6.1. Lineeirization of the two-dimensional stellar structure equations 

The dependent variables to be solved for are pressure Pt, temperature T, radius r, 
luminosity L (hereafter we use L to replace L', but remember that L is in solar units); the 
independent variables are chosen to be mass m (or s = Inm) and angular coordinate 6. The 
magnetic field variables Xj "^^j and 7 are also dependent variables. However, since we do not 
introduce their governing equations (such as the dynamo equations), we consider them to 
be given. All units arc in cgs, except for the luminosity that is in solar units. 

The construction of a two-dimensional stellar model begins by dividing the star into M 
mass shells and N angular zones. The mass shells are assigned a value Sj = logm^, where 
rrii is the interior mass at the midpoint of shell i. The angular zones are assigned a value 6j. 
A starting (or previous in evolutionary time) model is supplied with a run of (P/^-, T^j, r'^j, 
Lij, X^j, lij) for i=l to M and j=l to A^. 

Here wc take the general case B = {Br, Bg, B^) as the example to show how to solve 
the two-dimensional stellar structure equations. In order to write down the linearization 
equations, we introduce the following notations: 



6. 



Method of solution 



V 



Gw? p 



(187a) 



7^ 



m p 



(187b) 
(187c) 



r 



Anr^p pm ' 
PV, 
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C ^ ^^fe-T^'lA (187d 

Lq \ dt J Pm 



r = i3*V, 1,2,3, 10,---,15, 17 (187e 



2 de \ 2 d9 
^2 A cote dr'Y^ 

D" ^ (187j 



3Lq r'^Kpp. 



. ^ac mi , , , 

^ = or 2...- (1 + ^)^ (187P 



jr2 = 1 m pCpTlmV^^,,, (187q 

conv 

^3 ^ -^2^'^,, (187r 
Consequently, the stellar structure equations in the general case can be rewritten as follows: 

Qpf 15,17 



^ V+ ^'T^' + 0(2), (188a) 



ds 

1=1,2,3,10 
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ds 

9r' 
ds 
dL 
ds 



15,17 
1=1,2,3,10 



. e=i e=4 



39 09 d9 

where V = Vrad in the radiative zone, V = Vc in the convective zone. 



cot^ + 0(2), 



(188b) 

(188c) 
(188d) 



We calculate the derivatives of the dependent variables with respect to s by the central 
difference scheme, e.g.. 



but we simply use the difference scheme 



dP' 
ds 



+ 



dP' 

ds 



i-lj 



(189) 



r' — r' 

ij 



'89 



(190) 



to calculate the derivatives with respect to 9 because the first two of the two-dimensional 
stellar structure equations are singular at the poles. Thus, we can define a set of functions 
that should vanish at the solution of the stellar structure equations. 







PU) - 


F'^ 








= ir'- 


- 




= {Lij - 

2 


Li-ij) - 



1 

— i 

2 



15,17 



15,17 



€=1,2,3,10 



'i-lj; 



(191a) 

(191b) 
(191c) 



1 

—i 
2 



cot 9.: 



+E(4 + -^u))^Pi-Ji-. 



e=i 



+E(4 + ^ 



cot 9j 



''"a '"ii-i)}' 



e=4 



(191d) 
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where Asj = (sj — and i—2 to M, j=2 to A^. Vj^, • • • , D^?,^ are defined as follows 



pl7 



cot 6*0 , , 
Ir - — r 



ij-i) 



cot 6^ , , 
1 ^fr- 



■1 



1 - 



cot 6i 



2A^, ' 



cot 9i 



cot^ 

2A^/ ''^ 



Asi 



[(Xij - X'i-ij) + {p'ij - P'i-ij) + ilij - 



-^Kx'ij - x'ij-i) + {p'ij - p'ij-i) + {ilj - 

Mi - X^-y) + {p'ij - p'i-ij) + i^j - ^lij)], 

Mj - ^d-ij) + {p'ij - p'i-i,) + {^'Ij - ^'Uj)l 



ASi 
2C0t^ 

2Aej 



{{xiij - x:,-i) + (p^, - + {^"r, - 



We want then to solve for the set of (P/^, T^^ , r-^, L^ ) such that Fp — F!^ = F\ 
with x', 1?, and 7 specified. 



R 



(192a) 
(192b) 
(192c) 
(192d) 
(192e) 
(192f) 
(192g) 
(192h) 

(192i) 
(192j) 



The linearization of Eqs. (191a-191d) with respect to (5-P/;, ^T/^-, 5r[j, ^Lij) yields AMN— 
A(N — 1) — 4M equations for the AMN unknowns. The 2(A'" — 1) additional equations are 
supplied by the boundary conditions at the center. Prom Eqs. (172-173), we can define 



F'J = rl,-^[si-log(47rp„i,/3)], 



(193a) 



Fl^ = Uj-L^j- 



cot 9i 



£=1 



e=4: 



(193b) 



where j = 2 to A^. Another 2(A^ — 1) additional equations are supplied by the boundary 
conditions at the surface. From Eqs. (181a-181b), we can define 



(194a) 



Ff+^^ = L'j^^{\nLMj-a4P'Mj-a,T'j^j-ae), (194b) 

where j = 2 to A^. The 4M additional equations are supphed by the polar boundary 
conditions, 





= P' - P' 


(195a) 




rjif rjil 

= -'il ~ -'i2) 


(195b) 






(195c) 




= Lil — Lj2, 


(195d) 



where i = 1 to M. The F equations are linearized, 

where w — P,T, R, L; i — 1 to M; and j = 1 to A^. The summation over 1 has non-zero 
terms only for 1 = i-1, i; the summation over k has non-zero terms only for k = j-1, j. See 
appendix A for the matrix coefficients. 



6.2. Solution of the linearized equations 

Rather than solving the (4MA^)^ system of equations directly, we takes advantage of 
the specific form of the equations and especially of the large number of zero elements in the 
matrix. Prom Fig. 2 we can see that only 12 by AMN elements are non-zero at most. 

The matrix is reduced in a forward direction (i = 2 to M) as the coefficients are defined 
and is then solved in the backward direction (i = M to 1) for the corrections (^-P/,, 
SRlp SLij) for j = 2 to N. The reduction procedure begins: (i) For j = 2, to use the polar 
conditions at ^ = to eliminate those elements with subscriptions l = i, k = j- l (i.e.. Block 
III defined in Appendix A), which can be done by simply adding Block III to Block II. At 
the end of this step, the matrix equation for a specified j looks like Fig. 3a for a 4-point 
star in the mass coordinate (including the center and surface boundaries); (ii) to use the 
central boundary conditions to eliminate the first two columns in Block I for i = 2; (iii) to 
continue diagonalizing the four bottom rows for i = 2; (iv) to store the right-hand side and 
the elements in the rightmost columns, see Fig. 3b. After this reduction is completed, the 
bottom two rows of the first part of the coefficient matrix become the "central boundary 
equations" for the F equations of the next pair of mass points. The method is repeatedly 
applied until the surface is reached, whereupon the surface boundary conditions complete 
the set of 4M equations, see Fig. 3b-c. For the back solution (i) the values of ((^-Pm2> ^'^{12) 
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are first calculated, (ii) then the values of {5R^2-, SLi2, (^-P/_i2) <^T^-i2) for i = M to 2 are 
calculated using the stored elements of the array and {SPI2, STI2), (in) and finally the values 
of {6R[2, SL 12) are computed from the central boundary conditions and the values of {SP[2, 
STI2), see Figure. 3d-f. Since the sub-matrix with j = 2 has been diagonalized, we can use 
it to diagonalize the sub-matrix with j = 1 and 3. For j = 3, we use j = 2 as the "polar 
boundary conditions", and so forth. Finally, we solve the matrix equation, the results are 
stored in the right column. 

6.3. Advancing the model 

These routines arc based on the work of Prather (1976) and their revised implementa- 
tions in YREC (Pinsonneault 1988; Guenther et al. 1992; Guenther & Demarque 1997). 

6.4. Time steps 

In this section we use the cgs unit for luminosity L (erg/ s) and use X (Y) to represent 
the mass fraction of hydrogen (helium). The angular zone index (i.e., j) is 2. 

The timing routine calculates the time steps based upon a hydrogen- or helium-burning 
source. Let Lh (erg/s) be the total hydrogen-burning luminosity, and Lne, the helium 
luminosity. There are two time steps, 

A^H = hydrogen-burning time step, 
Att = total time step (i.e., for entropy and helium), 

where Att 7^ Atn only if the hydrogen shell is being shifted outward. If Lh = 0, the following 
section for hydrogen burning is skipped. 

For hydrogen-core burning {Xcoie > X^^^): a time step corresponding to a set of reduc- 
tion in Xcore is calculated. Let i be the innermost point if the core is radiative (i = 1) or the 
outermost convective point if the core is convective. Then, the change in Xcore is computed, 

AXeore = miujAX^^, A/r^ • X'} , 

and the time step is 

Ati = AtH = 6 • 10^« • AXeore " m'/ L\ 

where is the mass of the core (gm) and U is the luminosity of the core (erg/s, assumed 
to be mainly hydrogen burning). 
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When the core-burning criterion no longer apphes (Xcore < ^T^'^): a hmit is placed on 
the total amount of mass that may be burned, 

Am = A/m ■ Mq ■ Xenv, 
AtS = 6 ■ 10'^ ■ Am/La. 

If there is a hydrogen-burning shell {Xcore = 0), the timing routine locates it. Let the 
subscript o denote the inner edge of the shell (first point where X > 0); let 1/2 denote the 
mid-point of the shell (X = |Xenv); and let 1 denote the end of the shell {U — U~^ < 10~^-L 
or X = Xenv or en = 0). There is a limitation set on the maximum depletion at the mid-point 
of the shell, 

At](^ = AXi/2. 

With the exception of the core-burning phase, the new hydrogen burning time step is limited 
by the previous total time step, 

Atui new) = mm{ 1.5 • Ait (old), Ai^, At^^}. 

If there is to be no shell shifting then one sets Att — Atu- If the hydrogen shell is to be 
shifted outward through Am^ in mass, then the shift time step is computed as 

Aishift = 6 • 10^^ • Xi • Arus/Ln, 

and the total time step is 

Att = A^H + Ai,hift. 

If there is a hydrogen shell (Xcorc = 0), the helium burning is examined. For helium-core 
burning (Fcore > X™™ and L^ore > 0.1^^©), the maximum helium depletion is 

Ayeore = min{Ay--, A/-- . yeore}, 

and the helium time step is 

A^He = 5.85 • 10^^ • AFeore ' Mq/L^,^. 

For helium-shell burning (i^ore < X"^^^), the amount of mass burned through by the helium 
shell is hmited, 

A^He = 5.85 • 10^^ • A/,„ • Me/Lne- 
The helium time step places an upper limit on the previously computed hydrogen time step, 

Att = min{Ait, A^He}, 
A^H = min{Ait, A^h}- 
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The following parameters used in the determination of the time step are read in at the 
beginning of each model run. Their typical values are given as follows: 

X^^ = 0.001, AX^ = 0.04, AY^ra^'' = 0.02, 
A/-- = 0.5, A/-- = 0.3, A/^ = O.OOISMq, 

AX^^ = 0.10, Am, = 5 • IQ-^Mq. 

Of course, we can also use a fixed time step to advance the model. 



6.5. Composition advance 

The mixing routine performs all the operations on the model that are needed by the 
application of the time step to increase the age of the model. The routine first checks that 
there is no mixing within the hydrogen shell if the shell is supposed to be shifted. If there is 
such mixing, the shifting is suppressed (i.e., set Att — Atu). 

Each mass element (the mass contained in the zone defined by m e [mj_i,mj] and 
9 e [dj-i, dj]) is burned individually by computing the energy generation rates for the physical 
conditions existing in that mass element from the previously converged model. Since the 
program stores only the values of hydrogen, total metal and oxygen abundance, the change 
in these quantities is computed as 

X(ncw) = X(old) - (dX/dt) ■ At, 
Z(new) = Z(new) + (dY/dt) ■ At, 
Xi6(new) = Xieiold) - (dXo/dt) ■ At, 

where At — Att inside the hydrogen shell {X — 0) and At — Atn elsewhere. 

The routine then mixes those zones that it is instructed to do by being given a set of 
indices (i = ii to 12 and j = ji to J2), 

(k=i2 1=32 \ /k=i2 l=j2 

Y.^^'ki ■ XkA ■ lY^^akij . 
k=ii l=ji / \k=ii l=ji / 

The weights Uki are proportional to the amount of mass associated with zone kl are set up 
in the point readjustment routine. 

If the hydrogen shell is to be shifted, the routine calculates 



Asshift ^{S- Sy2 + - 5V4)/ In 10, 
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where 5 = Ams/mi/2 << 1. The points in the hydrogen shell are shifted by Asghift, 

So < Si < s i ^ Si (new) = Si(old) + Asshift, 

where Sj = logmj. The points up to a distance fs ■ Asghift in front of the shell are squeezed 
together, 

Si < Si < Send Si(new) = Si(old) + Send - Si{o\d)]/ fs, 
where Send = Si + fs- ASsMft- 

For all of these shifted and squeezed points the changes in P' and T' must be preserved 
for the calculation of the entropy energy term in the subsequent model. Thus for every 
Si(ncw), one must locate s;(old) such that S;(old) < Sj(ncw) < Si+i(old) and then interpolate 
linearly in s to get the old values of P' and T' which correspond to the new value of s. Then 
the effective changes are stored, 

AP^j = i^^(new s) -f';.(prc-shift s), 
ATI- = 7;;.(news)-7;;.(pre-shifts). 

For the region in front of the shell that is squeezed, it is desirable to preserve the original 
composition gradient if such a gradient exists. The values of X, Z and Xig are interpolated 
linearly in s as arc P' and T'. Note that the shifting process affects only the value of s and 
not the values of (P, T, R, L, X, Z, Xig) with the exception of (X, Z, Xig) in the squeezed 
region. 

The mixing routine finally checks on the physical sense of the new composition at all of 
the points, 

X = max{X,0}, 

Z = mm{Z,l- X}, 

Xie = max{Xi6, 0.99 • 10-^ • ZcNo}- 

The first two requirements are obvious; the third requirement brings the value of Xiq up to 
the approximate equilibrium value while turning off the Xiq burning rate that is calculated if 
XiQ > 10~^ • ^CNO- The value of ^cno — Z — Z^ where Z^ is the original weight abundance 
of all non-CNO metals. This method allows for the enrichment of CNO elements from the 
helium burning. 

6.6. Mixing zones 

Consecutive mass shells, which are determined to be convective (Vrad > Vad) in the 
previously converged model, are mixed together. 
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If there is a helium-burning convective zone, the semi-convective instabihty is treated as 
an over-shooting (Castellani et al 1971). The composition is first burned and mixed according 
to the standard convection zones. At the first radiative point outside a hehum convective 
zone, the quantity / = V™^/V^ad defined where the radiative gradient is computed with 
the (s, P, T, r, L) values of the radiative point and with the composition of both the radiative 
point {^^^) and the interior convective zone (™*). The original convective zone is extended 
outward through the radiative region for all the points at which /Vrad > Vad- 

This over-shooting region is restricted to the helium core (X = 0) and is limited by the 
condition of Castellani et al (1971) that defines a maximum radius i?max of the over-shooting 
mixing. 




(1 - l^{r)/i/ndr < (1 - V^VV;"4)L,At/(407rPei?c'), 



where the subscript c refers to the (s, P, T, r, L) values at the edge of the original convective 
zone. Here // is the mean molecular weight. The composition is then re- mixed from the 
beginning of the convective zone to the maximum extent of the over-shoot region. 



6.7. Point readjustment 

The point readjustment routine reflects all of the points between successive models. 
This routine starts with the central point and places each subsequent new point i so that all 
of the following criteria are met: 

P' — P' < AP' 

Li2 — Li_i2 < A/l • Lm2- 

All of the new values are interpolated linearly in s by locating the old point 1 such that 

Sz(old) < Sj(new) < s/+i(old). The fundamental variables (s, Pt, T, R, L), the composition 
(X, Z, Xig), the density and the entropy terms (AP', AT') are relocated between the center 
and outermost points for all angular zones. These variables are stored in temporary arrays 
and arc transferred to the original arrays once the process is completed. In addition to the 
and A/P^^ points remaining fixed, other points may be retained: 

1) the first radiative point (outer edge of convective zone), 

2) the innermost point of the convective envelope, 

3) the edge of the hehum core (X = 0), 
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4) composition discontinuities, X12 — ^^-12 > AXdisc or Z12 — ^;-i2 > AZdisc- 

The point routine then recalculates the weights assigned to each mass shell based upon 
the mass values at the preceding and following mid-points, 

rui = 10'% 

ai = ^(mi+m2), 

«i = ^("^j+i - m-i) for i = 2 to M - 1, 

qm = Mtot - ^{friM + friM-i)- 

The value rrii defines the location of the i*'^ shell, and is the number of grams contained 
in the shell. 

Additionally, the point routine adjusts the temperature of the outermost M*'^ point by 
adding a new point or deleting some old points. Given the desired temperature range Tmi„ 
to Tmax, then if Tm < T^in the outermost point j < M such that Ti > T = ^{Tyam + T^max) 
is selected as the new surface point. The points Z + 1 to M are deleted. If Tm2 > Tlnax the 
process is more complicated. The last atmosphere that was integrated will have stored the 
values of (Satm, -Patmj, T^tmj, raXmj) for the first inward integration step in which Tatmj > T. 
The new point M + 1 is added with the following values. 



-sm+i - 


~ "^atm 


P' - 

M+lj — 


P' 
atmj 


T' - 


T' . 

atmj 




/ 

' atmj 


Lu+ij = 




-'^M+lj = 


-- Xuj 




= Zuj 







6.8. Model calculation sequence 

The following list describes the sequence of calculations that is used in computing a 
series of stellar models. This sequence is the same for both one- and two-dimensional model 
calculations. 

(0) Input a model and compute a time step. 

(1) Locate the mixing zones and advance the composition and hydrogen shell for the given 

time step. 
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(2) Calculate element diffusion for the given time step. 

(3) Readjust the points in the mass coordinate in the model. This step is the main source 

of numerical errors and should be switched off for high precision calculations such as 
solar variability applications. 

(4) Calculate the entropy terms (AP' and AT'). Just zero them at the beginning, and give 

an estimate using their temporal change rate times the given time step. 

(5) Add the predictable corrections to (P', T', r', L) if their temporal change rates are 

available (after advancing one time step). This allows us to use a much larger time 
step and save a lot of computation time. 

(6) Specify the magnetic field configuration by selecting the functions x(m, 9), '&{m, 6) and 

7(m,e). 

(7) Retaining the old surface (or envelope) triangle and surface boundary conditions, do 2 

iterations for corrections to the dependent variables (P', T', r', L) and apply a given 
fraction (< 100%) of the corrections. 

(8) If necessary, relocate the surface triangle for the partially converged model and compute 

new atmospheres and surface boundary conditions. 

(9) Iterate until the model converges. 

(10) Refine composition and iterate until the model converges for solar applications that 
need a high precision. 

(11) Repeat (9) once for solar applications. 

(12a) If the corrections are excessively large at any time or if the model does not converge 
after many iterations (say 20), then retain the previous model that has been stored on 
disk and stop. 

(12b) If the model has converged, 

(i) compute a new time step, 

(ii) perform the requested printing, 

(iii) store the model temporarily on disk, overwrite the previous model, 

(iv) return to step (1). 
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7. Test 1: Two-dimensional standard solar model 

In this test, we investigate how different resolutions and different boundary conditions 
affect the two-dimensional solar models in the standard case (zero-magnetic field). 

Starting from a one- dimensional ZAMS (Zero Age Main Sequence) model, we move the 
fit point to the surface where the mass coordinate s = 1 ■ 10^^"^ from the usual location 
s = 1 ■ 10^''"' in a stair stepping way. The (ZAMS and the advanced) models are determined 
by the following parameters: the minimum and maximum change in s between Henyey grid 
points, 1 • 10~^^ < As < 8 • 10~^; the maximum change in w' {=P', T', r' and L/Lq) 
between Henyey grid points, is \5w'\ < 5.2834- 10~^. The convergence criteria for the stellar 
parameters are \SP'\ < 6 • 10"^ \ST'\ < 4.5 • 10"^ \Sr'\ < 3 • 10"^ and \S{L/Lq) \ < 9 • 10"^. 
The convergence tolerance on the rhs of the P and r equations is 3- 10~^, and the convergence 
tolerance on the rhs of the L and T equations is 2.5 ■ 10~^. We also require \SP'/P'\ < 9, 
\6T'/T'\ < 5, \6r'/r'\ < 5, and \6L/L\ < 90. All these requirements must be satisfied 
simultaneously when we apply the correction to the model. This is why we have to move the 
fit point in a stair stepping way. Otherwise, the correction is too large and the solution will 
diverge. The model has about 2401 grid points in the mass coordinate s, i.e., M = 2401. 
We also test the cases with M = 1201, 601, 301. 

When this one-dimensional convergence has been obtained, the angular part of the two- 
dimensional grid is selected. Unlike the mass coordinate s, which is not uniform, we simply 
equally divide the angular coordinate 9 in the range 9 e [0, 7t/2], 9j = {n/2){j — 1)/{N — 1), 
where j — 1 to N. We use the converged one-dimensional model for every angular zone. We 
use N — 10, 19, and 37 in this test. 

The solar mass is Mq = 1.9891 • 10^^ gm. The initial metal mass fraction is assumed to 
be Z = 0.022 at ZAMS. The model will evolve from ZAMS to the current age of the Sun (4.55 
Gyr). The hydrogen mass fraction and mixing length parameter (ratio of the mixing length 
over the pressure scale height) are determined by the requirement that the solar model 
at present reproduce the observed radius {Rq — 6.9598 • lO^'^ cm) and luminosity (Lq = 
3.8515 ■ 10^^ erg/s). We first use the one-dimensional code to generate a one-dimensional 
standard solar model as the reference. We then use the two-dimensional code to generate 
the two-dimensional zero-magnetic field solar models with different M and N and different 
surface boundary conditions and compare them with the reference. Our aim is to investigate 
if we can get a two-dimensional high-precision solar model. 
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7.1. Convergence 

First of all, convergence is the most important requirement in model calculations. There 
is an intrinsic divergence at the poles in Eqs. (188a-188b), and (188d), which results from 
the terms that contain cot 9. In order to solve this intrinsic divergence problem, we require 
both Eq. (185) and (186) at and near the poles. In practice, we zero Eqs. (195a-195d), where 
subscript '1' indicates the pole {9 — 0) and subscript '2' means the adjacent point to the 
pole. The denser the grid in the second dimension, the more severe the intrinsic divergence 
problem. Therefore, it is desirable to use fewer grid points in the second dimension for the 
sake of convergence. 

Since we have neglected the second-order derivatives with respect to 9 that are beheved 
to be smaller corrections to Eqs. (188a-188b), and (188d) than the first-order derivatives 
with respect to 9, we neglect those second-order derivatives to remove the divergence due to 
the numerical errors caused by them. 

There is a numerical divergence problem due to the possible equality between rjj, 
Pij, and Ty-i, Pij-i, respectively. When, say, rij equals Vij-i, the difference between 

them, TZ' = r^j — Vij-i, vanishes. In this case, the derivative of the difference with respect to 
(dTZ' / dvij) or rij^i {dlZ' /drij^i) should also vanish (i.e., dVJ jdvij = 0, or dlZ' / drij^i = 0, 
when Tij = Tij^i). If one sets dTZ'/dvij = 1 and dTZ' / drij^i = —1 no matter whether Vij 
equals rjj-i or not, one will run into a numerical divergence problem. We introduce the 6r, 
Sp, and St functions in the Appendix A to solve this divergence problem. 

The fourth divergence problem is due to the numerical error caused by numerical inte- 
gration of pm that affects the ratio p/pm, which is a two-dimensional correction factor that 
appears in all the stellar structure equations, Eqs. (188a-188d), noticing that the intrinsic 
singularity requires that the fewer the grid points at 9 the better. The numerical integral is 
usually made in terms of the trapezoidal rule, which is of the second-order in accuracy. De- 
upree (1990) adds more grid points to increase the integration precision when the numerical 
integral is performed. We find that it is more efficient to introduce a normalization factor in 
the integral, as shown in Appendix A. 2. 

When the radiative diffusion approxiamation (i.e., A = 0) is used, the code converges 
very well. This approxiamtion is not valid near the surface. If we use the temperature 
gradient at the surface to replace the actual gradient V^, the code also converges well. 
However, if we use the exact expression given in Eq. (55), we cannot get a converged model. 
The main cause is due to the numerical errors in the numerical derivatives associated with 
Ao. 
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7.2. Resolution 

If the convergence solves the internal- or self-consistency problem, then model resolution 
will address the external-consistency issue. Our reference model, i.e., the one-dimensional 
standard solar model, is almost the same as the best model described by Winnick et al. 
(2002), who emphasize its comparison with various observations. 

From numerical experiments using different resolutions in both dimensions, we find that 
the model is not sensitive to the resolution in the angular coordinate, but very sensitive to 
the mass coordinate, see Fig. 4. This figure compares 4 mass resolutions, in which the lower 
resolution is obtained by taking out one mass point every two points from the adjacent 
higher resolution model. Fig. 5 zooms in to compare the models with the highest and second 
highest resolutions. 

We compare different angular zones in Fig. 6 to make sure that the two-dimensional 
model is self-consistent in the angular direction. Fig. 7 shows that the two- and one- 
dimensional solar models with the same mass resolution are in very good agreement. 

7.3. Surface boundary conditions 

Until now we have used only the standard surface boundary condition used in YREC 
(Pinsonneault 1988; Guenther et al. 1992). If we use these standard model surface values 
of pressure and temperature as Deupree's reference values, as indicated in section 5.2.2, we 
obtain the same results, as seen in Fig. 8. The sohd lines use Eqs. (184a-184b). In order to 
investigate how errors in the reference pressure and temperature affect the result, we add 1% 
to Pfcf given in Eq. (184a) and 0.1% to T^f given in Eq. (184b). The result is shown by the 
dotted lines in Fig. 8. From the dotted lines we can see that errors in the surface boundary 
condition have larger influence on the outer layer than on the deep part of the model. 

It is inevitable to introduce some errors when Pref and T^ei are selected in model cal- 
culations. Nevertheless, Deupree did not need to worry much about it, since his interest 
focused on the core convection. In contrast, we should be cautious to use Deupree's surface 
boundary condition, because we want to apply to solar variability that takes place in the 
convective envelope. 

The model is less sensitive to the error in the reference pressure than to that in the 
reference temperature. 
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8. Test 2: Shell-like magnetic fields 

Shell-like magnetic fields depend upon only the radial coordinate r. Any physical mag- 
netic field should be free of divergence. The following magnetic fields are both radius- 
dependent and divergence-free: 

B = (0,0, /(r)), 
B = (C/r2,0,0), 
B = (C/r^0,/(r)), 

where /(r) is an arbitrary function of r, and C is an arbitrary constant. If we assume that 
there is no magnetic field in the radiative zone of the Sun, we have C = 0. Consequently, 
the unique physical shell-like magnetic field is 

B=(0,0,/(r)). (197) 

This is the case described in section 4.2, in which 

1 o cotedr'^ 



2 'V 2 de 

m. ^ 

B 



j4 ^ XP 



Comparing the 2-D stellar structure equations (Eqs. 124a-124d) with their 1-D counterparts 
(e.g., Li et al. 2003), we can see that the terms and/or factors in the brace are due to 2-D 
effects. 

In the solar variability applications, we use a standard solar model at the current age 
(t=4.55 Cyr) as the initial model. We apply a cyclic magnetic field to the model and use 
one year as the time step to advance the model. 

As in the one-dimensional case, we specify x ^ functions of time t (or sunspot number 
Rz) and the mass depth niD — log(l — m/Mo) as 

xifUD, Rz) = Xo(-Rz) exp[-|(mi:) - mT,cY/a\ (198) 

where rriDc specifies the location and a specifies its width, xo is determined by 

Xo{Rz) = ^{140 + [1 + log(l + Rz)fY, (199) 

where Bq is an adjustable parameter (unit: gauss), pc is the density at the mass depth of 
moc- In this case the magnetic variable-related derivatives reduce to 
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^ dlnx 9 In mo dlnm 

^ 5 In mo c)lnm dhiPT 

~ (t2 In 10 10™D GM|(1-10™D)2- 

In this test, x does not depend upon the angular coordinate 9, as required by a sheU-hkc field. 
The resultant models should be the same as we obtained in the one-dimensional counterparts 
(Li et al. 2003). The method of solution used in this study guarantees this test, as confirmed 
by actual model calculations. 



9. Conclusions 

A high-precision two-dimensional framework for treating stellar evolution with magnetic 
fields has been developed and successfully tested. The required high precision is achieved by 

(1) using the mass coordinate to replace the radial coordinate, 

(2) including the convection instability, 

(3) including a stellar atmosphere, 

(4) allowing element diffusion, 

(5) using fixed and adjustable time steps, 

(6) adjusting grid points. 

The code has the potential to include rotation and turbulence, but does not have the 
potential to generate them like a fully hydrodynamic code. 

We thank R.G. Deupree for many discussions and constructive suggestions during his 
stay in the department. We also want to thank Dr. Christian Straka for helpful discussions 
on many aspects of this paper. We wish to acknoledgc the anonymous referee, whose efforts 
have resulted in considerable improvements to the paper. This work was supported in part 
by NSF grants ATM 0206130 and ATM 0348837 to SB. SS and PD were supported in part 
by NASA grant NAG5-13299. P. V. was supported by Regione Lazio funds. 
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A. Coefficient matrix 

Eq. (196) consists of a set of non-homogeneous linear algebraic equations. If we use A 
to represent the coefficient matrix, and use B to represent the non-homogeneous term, this 
equation can be written down as follows: 



A - 5w ^ B, 



(Al) 



where 



Sw 



( Sr' \ 



11 
5Ti, 



6r', 
6L' 



Ml 



Ml 

8P'mi 
^T'mi 



(> I IN 

5L[n 



Sr' 
5L 



MN 



MN 



(A2) 
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is a column matrix. B is also a column matrix, 



V 



Ft 
Ft 



Ft 
Ft 
Ft 
Ft 



pM2 
R 

77M+I2 



Ff 

Ff 
Ff 



pMN 
R 

pMN 
jpM+lN 
pM+lN 



(A3) 



The coefficient matrix A has elements Only those elements with I — i — 1, i and 

k — j — 1, j may be nonzero, as shown in Fig. 2. We work out these nonzero elements in 
this appendix. 
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A.l. Useful partial derivatives 



The partial derivatives of the differential equations are required for the linearization. By 
defining the shorthand notation dxY = dY/dlogX, we can calculate the useful derivatives 
as follows. 

The following derivatives are almost the same as in the one- dimensional case (see Prather 
1976) except for those terms due to p/pm- These derivatives are nonzero for 1 = i-1, i and k 
= j. If k = j-1, they vanish. 



dpV 

dpTZ 

dpC 
dpC 



-A-V 
drV = 
-V 
-3•7^ 




= — CKr 



• 7^ 



5m ■ Tl 



Mr 



de 



Mr 



dlnP 


de 
dlnT 



+ 



+ 



dS 
dhTP 



dS 
dh[T 



/At 



/At 



In the convective zone, we have 

dR% 
dL% 
dp% 
dT% 

In the radiative zone, we have 

dpTr 
dL% 

dpTr 

dT% 



(ainVc/ainr-4) -T; 


(ainVc/91nPT-l)-re 

(9Ve/9inr)-r, 

Tr/L 

{d\nK/d\nPT)T ■% 
[{d\nK/d\nT)p-A]-Tr 



P 

Pm 

_P_ 

Pm 



The formulas for the various partial derivatives of the physical quantities will be presented in 
the following subsections. The equation of state calculates p, a, 5, Cp, Vad and the pressure 
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and temperature derivatives of these quantities (see section B.l). Energy generation rate e 
is a function of p and T, too. So {de/dT)p and {de/dP)T can be expressed by {de/d\nT)p 
and de/ din p)t (see section B.2): 

{de/d\nT)p = {de/d\nT)p + {de/d\n p)T{d\n p/d\nT)p 
{de/ din Pt)t = {de/dln p)T{dln p/dlnT)p 

The derivatives of the convective gradient Vc which are presented in section B.4. 

The entropy term contains the only exphcit reference to any time-dependence in the 
stellar structure equations. It can be reformulated as follows: 

S = -(PT5/p)-(AT7Vad-AP') 
{dS/d\nT)p = S[6 + {d In 6 /d\nT)p]-{P6/pV^)[l-{dlnV^/dlnT)p- AT'] 
{dS/d\nPT)T = S[l-a + {d\n6/d\nPT)T\ 

+ (PVp)[1 + (ainVad/a In Pt)t • AT'/Vad] 

where (AP', AT') are the changes between successive models. 

The following derivatives are new. Similarly, these derivatives are nonzero for 1 = i-1, i 
and k = j. When £ = 1, we have 

dpB^ = -A-B"" 

dLB' = drB'^O 

dpB^ = -B' 

dR%^ = • (91nVe/91nr-4) 

9lT; = 

dpT^^ = T^^ -[d InV c/d In Pt- I] 

drrT^ = ■{dXnVjdlnT) 

dnT^ = -4-r/ 

diJr = ^V^ 

dpT^ = -{d In K/d In Pt)t 

dT%^ = r/ • [(a In In T)p^ -4] 

When £ — 2,we have 

dpB^ = ~B^ 

dLB^ = drB^^Q 
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When £ — 3, we have 



dpB^ 
dpT^ 

da: 

drX' 



hi 

drB^ 



dpB^ 



dpT! 

da! 

dpT! 

da! 



dpT; 



3 _ 

da! = 

dpT! = 

da! = 

When e = 10, 11, 13, 14, we have 

dpB' : 
dpB^ - 
drB^ - 

dpT! 
da! 

dpT! 

da! -- 



-B' 

-{dlnWc/dlnR-l) 



-[dlnWc/dlnP -1] 
■ (91nVc/91nr) 

-r! 
r!lL 

T! ■ (a In K/a In Pt)t 

T! ■ [{dhiK/d\nT)p^-A] 



-3-B' 


-B^ ■ am 

B'-6m 

■ {d\nVc/dlnR-3) 



-(dlnVc/dlnP-am) 

■ (dlnVc/dlnT + 5m) 

-3 • t! 

V/L 

T! ■{d{\YiK/d\TlPT)T-am + l) 
T! ■{dhlK/dhlT)p^ + 5m-4) 



dLB' 
7i I 



13 

-B'-6 



= 

a-1) 



T! -{dhiVjdXur) 


T! ■{d\uVjd\nP + a-l) 
T! -{dluVjdluT -5) 
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OtT! 



drTS 



When £ = 12, 15, 17, we have 





Tl ■[{d\nK/d\nPT)T + Oi] 

t;^- [(ain«/9inr)p^ -5-4] 



dpB^ 

da! 

da! 



da! 



dpT! 

da! 



-3-B' 
OtB^ = 
-B' 

T! ■ {dlnVjdlnr -3) 


T! -{dlnVc/dlnPT-l) 
T! ■ (dlnVc/dlnT) 

-3 ■ t! 

r!/L 

T! ■{d\nK/d\nPT)T 

T! ■[{dlnK/dlnT)p^-4] 



When k = j-1, all derivatives of B's and T's vanish. 

Wc also need similar derivatives of Vs. When 1 = i-1, k=j, all derivatives of V, V, 
V, and V^ are zero. The nonzero derivatives are: 



10 



dp ,.V 
dp , .V^^ 



dp , = 



/As, 



-dp...V 



pl4 
10 



dPi-,, 
dr , V 



14 



-dr. 



When 1 = i, k=j, the nonzero derivatives are 



f cot 9j 

cot On 

H 

2Ae, 



cot 9i , . 
-{r- 



1 - 



cot 9j 
2Kf, 



-l) 



5r, 
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dR.V'' 



Or V 



dR.V^^ 
dR.V^' 

dp.v^^ 
dp.v^^ 



dp.V^^ 



14 



drV 



cot 6j 

f cot 6j 
1 



1 S^r/ _ / ^ 



(-^7 -^ii-l) 



1 - 



5r, 
2 AO, ^ 



1 



cot^ 



-dR^P', 



2A^^- ' 



cot 
2A^ 

-9^,P^'^(r^ - r^_J/A^„ 

P^a^/p/A^,-, 

P^ccy/As,, 

-V^aij/Asi, 
^9 cot ,■ 

2Aej ^ 

-V%j/Asi, 

-dT,V'%r[,-r[,_,)IAe,, 

-VdijSp/Aej, 

-V''5ij/Asi, 
V%j/Asi, 



-V 



2Ae, 



where 5r = 1 when rij — Vij-i ^ 0, and Sr ~ when Vij — rij-i = 0. 6p and 6p have the 
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similar meaning. When 1 = i, k=j-l, the nonzero derivatives are 
for £ = 1,2,3, 10, •••,15, 17, and 

dp.. ,V^^-dp.V^ 

for £ = 3, 12, 15, and 

-I — 1 -1 13 

for £ = 12, 15. 

We calculate the derivatives of and by taking the advantage of /mf^conv ~ const. 
The nonzero derivatives are listed as follows for / = i — 1, i and k = j. 

dnJ^' = -2-T\ 
dLT^ = 0, 

dpT^ = -JFi. [(91nK/91nPr)T + Q; + Q;J, 
d^T^ = :ri.[4-(91nK/91nr)p^ + (5 + 5^], 
BrT'' = -2•J^^ 
9lJ^' = 0, 

dpT'^ = (9 In Cp/a In Pt)t 

-/3[2q; + (ainK/ainPT)T + (91nCp/ainPT)T]}, 
d^T'^ = JF^. {1 + (a In Cp/a In T)p^ 

+/3[3 + 25- (a In In T)p^ - (5 In Cp/<9 In T)p^]}, 

BrT'' = -2•JF^ 

^lJ^^^ = 0, 

5p^3 ^ -ap^2.^,^^^3.^^1^y^^/^l^p^)^^ 

« x-4 _ x-1 / ^ + (91nVe/ainr-l)^^4 (convective), 
^^•^ - ^^-^ ■{ (radiative). 



•4 



(convective), 
I L (radiative). 



a ^ a I ^ + {d\nVjd\nPT + a-l)-T' (convective), 

a ^ c) I +(91nVe/91nr-5)^J^^ (convective), 

^ ^ ■» '''^^ + [(9 In «/9 In T)p -5-4] •JF^ (radiative). 



GmpV c. 
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= 0, 




= dpT^ 












= 0, 











^^^ + (91nVe/91nr-l)•^^ 



GmpVc ^ Ve/91nPT + a - 1) • ^^ 



Gmp _ ^ 
rPr ' 



rPr 
Gmp 



- 5 • Jp•^ 



rPr 

where /3 = (VconvAo)/(l + ^^conv/^^o)- 

A. 2. Numerical integrals 

The quantities p^, oim and bm are integrals over ^: 

\ \ r 

Pm{m,9) = — ■-/ dOr'^ {m, 9) p{m, 9) sin 9, 
^ ^ Jo 

dhipm\ Jq d9r'^ {m, 9)p{m, 9)a{m, 9) sm9 



a In Pt/t d9r^m, 9) p{m, 9) sin 9 ' 

51n \ Jo d9r'^{m, 9)p{m, 9)5{m, 9) sin 9 



am{m) = 

^ V^lnTy/p^ /;(i^r2(m,^)p(m,^)sine 

Of course, the luminosity L is an integral, too: 

L{m) = 2n d9r^{m, 9)Fr{m, 9)sin9^- / d9L'{m, 9)Lq sin 9, 
Jo 2 Jq 

where L' — Anr'^Fr/ Lq. In the one-dimensional case, we know the relationship on the solar 
surface: 

L = 47ri?VT4. (A4) 

If we define 



= l [ d9r\Mu,„9)sin9, 
^ Jo 

T^s ^ ^£d9r'T\M,,„9)sin9, 
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Eq. (A4) holds well in the two-dimensional case, where Mtot is the total mass of the star. 
We use the trapezoidal rule to compute these integrals. For example, 



p'^ = 

rm 






m- 


- Oi-i) 


L 




:2 W^'mi sin + ^Mt-x sin 




I — ^e-i) 






Y!i=2 \ (sin + sin 


9 — 6^. 


-i) ' 




1^1= 


a \ i^lie sin 0^ + r^<,_i sin Og. 




- Oe-i) 




Y!i^2 |(sin Ot + sin 6'^_i) (6'^ 




i) ' 



where is the total grid number in the second dimension 6. We have introduced the 
normalization factor [X]^2 |(sin + sin^£_i)(^^^ — to remove the discrete error. The 

other three integrals do not need the normalization factor: 

_ Yle=2 WitPi^^^it sin + r^ViPtf-ia^^-i sin g^-i) {0^ - O^-i) 

J2f=2 \ i'^uPie sin Oe + rf^^iPu-i sin 9e-i){0e - 9e-i) 
_ J2e=2 WuPii^ii sin + rf^_ipu^i5u-.i sin 9i-i){9e. - 9e-i) 
Y!!=2 WuPii sin 9e + r'^i_iPie-i sin 9e-i){9e - 9i-i) 
rp4 ^ E^2 K^mA sin 9i + rlje_iT^^_^ sin 9i-i){9i - 9e-i) 
T!!=2 sin 9i + rl^^_^ sin 9e-i) {9^ - 9i_i) 

because they have already had their own normalization factors. 

A. 3. Interior points 
A. 3.1. w^P 

There are three blocks in this group. They are: 

(I) l = i-l,k = j; 

(II) l = i,k = j; 

(III) 1 = i, k = j-1. 

We present the results one block by one block using the derivatives given above. 



3,12,15,17 

e=i 



AsidpVi-ij - 1 

15,17 



■As. 



i=l,2,3,10 ^=10,11,13,14 
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As,: 



.3,10,11,13,14 ^=10,11, 13,14 



^=10,11, 13,14 



--AsidnVij 



--Asi 
2 



'3,12,15,17 15,17 

Y ^4-^^+ E (Ci.+4)^.^^ 

. e=l ^=1,2,3,10 



--AsidpVij + 1 



--Asi 
2 



1 



15,17 



15 



Y dpBl,.V^+Y(^U + ^'.)9p^P' 

=1,2,3,10 ^=3,10 



-AsidrVij 



--Asi 
2 



15 



=10 



-67- 



Block III: 



_dFl_ 
dLij_i 

dP' 

_dFl_ 
dT' 



15,17 



^=1,2,3,10 



= 



^=3,12,15 



^=12,15 



A.3.2. w=T 
There are three blocks in this group, too. 



Block I: 











dF^' 




-li 


e-- 


dFi^ 
dPi-ij 


—AsidpTi. 


1 / 





dF^' 



i-lj 



15,17 



15,17 



£=1,2,3,10 



':ASi 



15,17 



i=l,2,3,10 £=10,11,13,14 



--AsidT%-ij - 1 

^ r 15,17 

--As. 



i=l,2,3,10 £=10,11,13,14 



Block II: 



dFl_ 
dR' 



1 



AsidR%. 



V 
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dFl_ 
'dL~ 



^3 



dFl_ 



dFl_ 
dT . 



=1,2,3,10 ^=1,2,3,10 
. . 15,17 



^=1,2,3,10 



--AsidpTij 



--Asi 
2 



1 



15,17 



15 



.€=1,2,3,10 €=3,10 



--As,dTT,j + 1 



2 ' 



15,17 



15 



.€=1,2,3,10 €=10 



Block III: 



_dFi_ 



=1,2,3,10 



= 



€=3,12,15 



€=12,15 



In this group only the first two blocks are nonzero. 



Block I: 



--AsidnTZi-ij - 1 



9L, 



= 



i-lj 
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dFl_ 
9TUj 



--AsidpTZi^ij 



--AsidTTZi-ij 



Block II: 



dFl 

dFl 
dLij 
dFl 

dPij 
dFl 



--AsidRTZij + l 



--AsidpTZij 



--AsidrTZij 



A.3.4. w ^ L 

Similarly, in this group only the first two blocks are nonzero. 



Block I: 



dFl^ 



dR' 



dFl^ 



dL. 



_dFl_ 
9PUj 



1 Asi cot 9j 
'2 A9-i 



.e=i 



1 As,- cot di 



-1 - 



2 A^,- 



--AsidpCi_ij 
1 As,- cot 9-i 



2 A9j 



e=4 
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dT! 



--AsidxCi^ij 



1 Asj cot 6j 

2 Mj 



1=1 



ij ij-l) 



£=4 



Block II: 



dFl 
dLij 
dFl 
dPij 



1 Asj cot 9j 

2 AOi 



1 As,- cot 9i 



= 1 - - 



2 A^,- 



--AsidpCij 

1 Asj cot 
~2 AOj 



--AsidxCij 

1 Asj cot 6'j 
~2 AOi 



where 5p — 1 when Pjj 7^ Pij-i-i ^iid 5p = when Pjj = Pij-i- The definition of 5t is 
similar. 



Block III: 



dFl_ 



1 Asj cot 9j 
2 
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dL 



ij-l 



_ lASiCOt^ 3 3 



where 5r = 1 when 7^ ^^y-i, and 5r = when = rij_i. 

A. 4. BoundEiry points 

A.4.i. Center: cu — r 
Central boundary points have only block II for w — r: 



= 1 







dLij 

dFH _ 1 
dP[^ ~ 3 

dT[. 3 ' 



mlj 



A. 4. 2. Center: u ^ L 
Central boundary points have blocks II and III for uj — L: 

Block II: 

2 



dFl^ _ _cotej 



1=4 
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dFl 

dFl_ 
dP' 



Block III: 



-dpCij 



cot 9j 



i=4 



= -drC 



cot 9j 





cot 9j 








= 


dLij_i 


dFl^ 


cot 6j 






dFl^ 


cot 6j 




AO, 



A. 4. 3. Surface 

Standard The surface boundary conditions are linearized as follows: 

5r'^^ + 0-5LMj-ai5P'j^^-a25T',^^ = -F^-"'', 
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Deupree's His surface boundary equations are simpler: 



where 



rpM+lj _ pi _ P' 

— ^M+lj -^ref' 

pM+lj _ rpl rpl 

— -'-M+lj -'ref- 



A.4.4. Pole 

The polar boundary equations arc extremely simple: 

sp:,-6p:, = 
st:,-st:, = 

5rl,-5rl, = 
5Lii - 5Li2 = 



B. Input physics 

B.l. The equations of state 

When a magnetic field is present, the equation of state relates the density p to the 
pressure P, temperature T, magnetic energy per unit mass x, the ratio of specific heats 7, 
and the chemical composition: 

P = p{.Pt,T,x;X, Z), 

where P = Pq + Pr + P„i is the total pressure, Pq the gas pressure, Pj. = the radiative 

pressure, Pm = XP the magnetic pressure, X the mass fraction of hydrogen, Z the mass 
fraction of elements heavier than hehum (the so-called metal mass fraction) . Its differential 
form is 



where 



dp (IPt ^dT dx 

— — ct— — ^ — , 

p P T X 



a — (dlnp/dlnP) at constant T, x, 
6 — —{dlnp/dlnT) at constant P, x, 
u — —{dlnp/dlnx) at constant P, T, 
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here X and Z are assumed to be constant. 

Since it is tedious to accurately calculate the equation of state from first principles, the 
equations of state are usually provided by the numerical tables as functions of (p, T, X, Z) for 
Po, -S" (entropy), U (internal energy), (dU/dp)T, Cy = {dU/dT)p, Xp = (d\nPo/dp)T, Xt = 
{d\nPo/dT)„ Ti = {dlnPo/dlnp)s, T'^ = r2/(l-r2) = 1/Vad, and r'3 = {d\nT/d\np)p,- 
1. The equation of state for the gas is taken from Rogers, Swenson, & Iglesias (1996). In 
order to take into account a magnetic field based on the EOS tables, one can use the following 
correction method: 

(i) Using the total pressure P — Po + Pr + Pm, the total internal energy U — Uo + 3Pr/ p + X: 

and the total entropy S — So + 4:Pr/pT + x/T to replace the gas pressure Pq, the gas 
internal energy Uo, and the gas entropy 5*0 respectively when interpolating to obtain 
the density for the given P and T; 

(ii) Using {Pq + Pm)/P to rescale Xp', 

(iii) Using Pq/P to rescale xt from the EOS tables and add iPr/P; 

(iv) Adding UP^/T to from the EOS tables; 

(v) Computing T'^ = Pxr/cvpT, Ti = Xp + Xt^'s: and r'2 = Ti/F^. 

Taking these as known, we can calculate 

a=l/Xp, S = Xt/Xp, i^ = Pm/P, Vad = l/r's, Cp = P6/pTV^. 
These quantities are used in calculating the convective gradient Vc. 



B.2. Energy generation 

The calculation of the energy generation includes the individual rates for the PP-chain 
(PPI, PPII, PPIII), the CNO-cycle with a simplified NO approach to equilibrium. The 
coefficients of all of the reaction rates and the formulae for most of them are taken from 
Fowler, Caughlan and Zimmerman (1975). 

The reaction rate for the PP-chain is actually that for the H^(p, e+z/)D^ reaction and 
assumes that all the other reactions in the chain are relatively instantaneous. The burning 
rate is then 

(dX/dt)pp = 4.181 • 10-^V^^r9"^^^exp(-3.380/T9/^) 

•0(q;') • (1.0 + 0.123^9/^ + l.OgTg/^ + 0.93871)) sec-^ 



-75- 



where Tg = ° K, the screening factor fg is set equal to 1, 

= l + a'[{l + 2/a'y/^ -1] 
a' = 1.93- 10^^(F/2X)2exp(-10.0/T9/^). 

The total energy of the PP-chain (subtracting the energy of the neutrinos which are pro- 
duced) is 

epp = 6.398 • 10^^'ilj{dX/dt)pp erg/gm/sec, 

where 

^ = 0.979/1 + 0.960/ii + 0.721/in 

/i = [(i + 2/«')'/'-l]/[(l + 2/a')'/' + 3] 

fn = (i-/i)/(i + r) 

/ill — 1 — /l — /ll 
r = 10^^-^'^^[X/(l + X)]T~'/'exp(-10.262/T9/'). 

The derivatives of epp can be found directly: 

((91nepp/(91np)r = epp 

(ainepp/ainT)p = epp[-2/3 + 1.1267/T9/^ + (91n0/(91nT)p + (ain^/ainT)p 

+ (0.041X9/-'^ + 0.727Tg/'^ + O.938T9) 
/(I + 0.123X9/^ + 1.09X9/^ + O.938T9)] 
{d\n(j)/dlnT)p = (2/0 -l)(l + 2/a')"i/23 333/^1/3 

{d\ni/j/d\nT)p = ^-^{[0.258 - 0.239/(1 + r)]((9/i/(91nT)p 

-o.239/in/(i + r)(ainr/ainr)^} 

{d\nT/d\nT)p = -1/6 + 3.4207/T9/^ 

{dfi/dlnT)p = -4{Q;'(l + 2/a')'/'[(g + 2/a')'/' + 3]T'-3.333/rJ/' 

In the calculation of the CNO bi-cycle, CN equilibrium is assumed and the CN cycle is 
assumed to be the only source of energy. The hydrogen-burning rate due to the CN cycle is 
then 

{dX/dt)cN = 1.202 • 10V^XNr9"^/^exp(-15.228/r9/^) sec'^ 
and the energy produced is 

5.977- 10^*(ciX/cii)cN erg/gm/sec. 
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The value of Xn (N^^ abundance by weight) assumes that all the carbon and nitrogen is in 
the form of N^^ 

— Z — — Xq, 

where Z is the total metal abundance by weight, Zm is the weight abundance of all non-CNO 
metals, and Xq is the weight abundance of O^^. The approach to NO equilibrium is taken 
as a simple burning rate of O^^ assuming O^^ equilibrium, 

dXo/dt = 9.54 • 10V^XoT9~^^/^^exp(-16.693/T9/^) 

-1.6 • lQ-^{dX/dt)cN sec"^ 

Between successive models the value of Xq is decreased at a rate of (dXo/dt) per second, 
and thus the value of Xn is correspondingly increased. Here are the derivatives of the CN 
energy production: 

((9ecN/<91np)T = ecN 

(aecN/ainT)p = ecN(-2/3 + 5.076/Tj/'. 

B.3. Radiative opacities 

An estimate of magnetic effects on the radiative opacities [k = k(T, p, X, Z)] can be 
found in Li and Sofia (2001). Since they are small, we use only the OPAL opacities tables 
(Iglesias and Rogers 1996) together with the low-temperature opacities from Alexander and 
Ferguson (1994). For X and Z the hnear interpolation is used, but for T and p the cubic 
splint interpolation is used. The cubic splint interpolation scheme allows one to obtain the 
derivatives of k with respect to T and p. These derivatives are needed in the linearization 
of the equations of energy transport. 

B.4. The convective gradient and its linearization 

The calculation of the convective temperature gradient Vc in the envelope of the stellar 
models employs the mixing length theory (Henyey et al. 1965; Lydon and Sofia 1995) when 
magnetic fields are taken into account. 

Defining 5' = Vrad — V'g^^, the Schwarzschild (1906) criterion is used to determine con- 
vection: S > means convective. In the deep interior convection zones Vc is set equal to the 
adiabatic gradient V^. 
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In the envelope the evaluation of Vc is more complex, and we solve 

F{x) = a3X^ + + oix - 1 = 

for y > such that F{y) = 0, where 

ai = V^as(f)iKTyCp){Hp/g55'y/^ 
3 

03 = -(pUj'^/V. 

We have defined 5' = (Vrad — V!^^), (f) = (1 + icj^)"\ u = pulm, and = 16y/2a. The root 
y is guaranteed to lie in the interval (0, +1) since F{0) = — 1 < and F{1) = ai + a2 > 0. 
Further, this root y is unique since the derivative of F, F'{x) = Sa^x^ + 2x + ai. is positive 
definite for x > 0. An initial estimate of the root y is made and a second-order Newton- 
Raphson correction is applied: 

Ay = -F{y)/F'{y) - ^[F(y)/F'(y)]^ • F"iy)/F'{y). 

The initial estimate of y is y = 1/ai, unless 03 > 10^ in which case y = (1/03)^/^ that follows 
the asymptotic behavior of the solution in either limit. Given the solution y, the convective 
gradient is computed: 

^c^V^ + {V,^-VJy{y + a,). 

The linearization of the convective gradient is cumbersome but can be calculated. We 
consider derivatives with respect to InP^, InT, Ini? and L. 



dhix dlnx dhix 

where 



(2y + «i)^i— + 



5 In X dhix 

85' _9V;,d 9Vad 



9 In X 9 In X 9 In X 
The derivatives of come from the equation of state, and are non-vanishing only for 
X — Pt or T. The derivatives of 

Vrad = (3/167rac)(KLL0PT)/(G'MT^) 

are non- vanishing for x = Pt, T, or L, 

«9Vrad „ f d\nK 

— Vrad I t + 



dlnPT V dlnPT 

(9Vrad _ / , (9 In K 



dlnPT ^"""^y ^^dlnTyp 
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In the radiative zone, the actual temperature gradient is equal to the radiative temperature 
gradient V = Vrad, its derivatives are given here. The opacity tables provide logK vs. (logp, 
log T) . In order to calculate 

{d\nK/d\nT)p = {d\n K/d\nT)p + {d\n n/dln p)T{d\n p/d\nT)p, 
((91n/t/(91n P)t = id\nK/d\np)Tid\np/d\nT)p, 

one needs {din K,/dlnT)p and {din n/dln p)t (see Iglesias & Rogers 1996, and Alexander & 
Ferguson 1994). The derivatives of y are functions of oi and 03, 

dy 1 / 91nai o da^, 
- ' aiy— h asy 



dlnx F'{y) \ dlnx d\nx 

which need the derivatives of Oi and 03, 

dlnai dlnd) dlnn dlnT dlnC^ IdlnHp 1 85' 
+ -TT. K3— — — ^ + - 



dlnx dlnx dlnx dlnx dlnx 2 dlnx 25' dlnx 

1 ding ldln5 

2 dlnx 2 dlnx 

dlna% dlnuj din 6 dlnai 

= 2 \ -. 

dlnx dlnx dlnx dlnx 

The derivative of 5 = —{dlnp/dT)p and Cp are computed by the equation of state. By 
calculating the derivative of 0, 

dlncf) 2 2^9lnuj 
dlnx 3 dlnx^ 

and by expressing Hp and g explicitly, one can obtain 

91nai 2 2^f^liia; dlnCp 1 d5 19 In (5 

dlnx 3 9 In X dlnx 25 dlnx 2 dlnx 

Pt dlnT dlnr dlnn Idlnp 
+ 3— + 2— — + ^ 



In a; dlnx dlnx dlnx 2 dlnx 
9 In 03 ^^dlnoj 91nai 

dlnx dlnx dlnx 

The derivatives of u with respect to Pt, T, r are straightforward, 

d In uj f dlnn 

1 + 



dlnPr \dlnPT/T 
d\nuj f dlnK 



dlnT \d\nT 
dlnuj 
dlnr 



p 



-79- 



REFERENCES 

Alexander, D. R., & Ferguson, J. W. 1994, ApJ, 437, 879 

Balmforth, N. J., Gough, D. O. & Merryfield, W. J. 1996, MNRAS, 278, 437 

Basu, S., Pinsonneault, M. H., & Bahcall, J. N., 2000, ApJ, 529, 1084 

Boothroyd, A. I., and Sackmann, I. J. 2003, ApJ, 583, 1004 

Callen, H. B., 1966, Thermodynamics (New York: John Wiley & Sons), p. 242 

Castellani, V., Giannone, P. & Renzini, A., 1971, Ap. & Space Sci., 10, 340 

Deupree, R. G. 1990, ApJ, 357, 175 

Endal, A.S., Sofia, S. & Twigg, L.W. 1985, ApJ, 290, 748 

Fowler, W. A., Caughlan, G. R., & Zimmerman, B. A., 1975, Ann. Rev. Astro. Astrophys., 
13, 112 

Guenther, D. B. & Demarque, R 1997, ApJ, 484, 937 

Guenther, D. B., Demarque, P., Kim, Y.-C. & Pinsonneault, M. H. 1992, ApJ, 387, 372 
Henyey, L. E., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 
Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 
Krishna Swamy, K.S. 1966, ApJ, 145, 174 

Li, L. H., Basu, S., Sofia, S., Demarque, P. & Guenther, D. B. 2003, ApJ, 591, 1267 

Li, L. H., Robinson, F. J., Demarque, P., Sofia, S. & Guenther, D. B. 2002, ApJ, 567, 1192 

Li, L. H. & Sofia, S. 2001, ApJ, 549, 1204 

Lydon, T.J. & Sofia, S. 1995, ApJS, 101, 357 

Pinsonneault, M. H. 1988, Ph.D. thesis, Yale university 

Prather, M. J., 1976, Ph.D. thesis, Yale University (Appendix A) 

Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902 

Sackmann, I. J. and Boothroyd, A. I. 2003, ApJ, 583, 1024 



-80- 



Schwarzchild, K. 1906, Gottinger Nachrichten, p. 41 

Sofia, S., Basu, S., Demarque, P., Li, L. H., & Thuillier, G. 2005, ApJ, 632, L147 

Sofia, S., Heaps, W., Twigg, L. W. 1994, ApJ, 427, 1048 

Sofia, S., O'Keefe, J., Lesli, J. R., and Endal, A. S. 1979, Science, 204, 1306 

Spruit, H. C. 1991, in Tiic Sun in Time, cds. C. Sonett, M. Giampapa and M. S. Matliews, 
University of Arizona Press, Tucson, p. 118 

Spruit, H. C. 2000, Space Sci. Rev., 94, 113 

Tlioul, A. A., Balicall, J. N., & Loeb, A., 1994, ApJ, 421, 828 

Turcotte, S., et al. 2001, in lAU Colloq. 185, Radial and Nonradial Pulsations as Probes 
of Stellar Physics, ed. C. Aerts, T. R. Bedding, & J. Christensen-Dalsgaard (San 
Francisco: ASP), 258 

Unno, W. & Spiegel, E. A., 1966, PAS J, 18, 85 

Winnick, R. A., Demarque, P., Basu, S., & Guenther 2002, ApJ, 576, 1075 



This preprint was prepared with the A AS IM^jX macros v5.2. 



-81 - 




5650 



5600 L 




_i I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I i_ 



2 5 

Age (Gigo year 



4 



5 



Fig. 1. — Reference values of pressure and temperature at the surface as functions of age. 
The dotted hues are polynomial fits to the calculated model (solid lines) using the standard 
surface boundary condition. 
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Fig. 2. — Linearization equation for a 4 x 4-point star. 
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Fig. 3. — Schematic Henyey solution for a 4-point star. The matrix block is denoted by 
O's, I's, and X's are non-zero. The right hand side is denoted by A, the elements changed 
through pivoting, by Y and B. The final reduction to the identity matrix is not shown. 
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Fig. 4. — Relative changes of pressure P, temperature T, radius R, luminosity L, density 
p, nuclear energy generation rate e, opacity k,, and sound speed C in two-dimensional solar 
models with different mass-coordinate resolutions (M = 2401, 1201, 601, 301) with respect 
to a one-dimensional standard solar model as functions of mass coordinate. 
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Fig. 5.— Same as Fig. 4, but only for M = 2401 and 1201. 
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Fig. 6. — Pressure P, temperature T, radius R, luminosity L, density p, nuclear energy 
generation rate e, opacity k, and sound speed C at different angular coordinates as functions 
of mass coordinate. 
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Fig. 7. — Pressure P, temperature T, radius R, luminosity L, density p, nuclear energy 
generation rate e, opacity k, and sound speed C in both one- and two-dimensional solar 
models that have the same mass resolution as functions of mass coordinate. 
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Fig. 8. — Relative changes of pressure P, temperature T, radius R, luminosity L, density p, 
nuclear energy generation rate e, opacity n, and sound speed C in the two-dimensional solar 
model with Deupree's surface boundary conditions [sohd line: Eqs. (184a-184b); dotted line: 
1.01-Eq. (184a), 1.001-Eq. (184b)] with respect to a one-dimensional standard solar model 
as functions of mass coordinate. 



